提交scripts

libs/Math/TRWS/CHANGES.TXT

@@ -0,0 +1,13 @@
+Changes from version 1.2:
+
+- Fixed bug in typeBinaryFast
+- Added the option of reading min_marginals
+- Added function AddRandomMessages()
+
+Changes from version 1.1:
+
+- Fixed bug in memory allocation
+
+Changes from version 1.0:
+
+- Modified syntax of 'friend' declarations to make it compile under unix

libs/Math/TRWS/LICENSE.TXT

@@ -0,0 +1,31 @@
+Tree-reweighted max-product message passing algorithm (TRW-S), verstion 1.3.
+
+--------------------------------------------------------------------------------
+
+This Microsoft Research Shared Source license agreement ("MSR-SSLA") is a legal agreement between you and Microsoft Corporation ("Microsoft" or "we") for the software or data identified above, which may include source code, and any associated materials, text or speech files, associated media and "online" or electronic documentation and any updates we provide in our discretion (together, the "Software"). 
+
+By installing, copying, or otherwise using this Software, found at http://research.microsoft.com/downloads, you agree to be bound by the terms of this MSR-SSLA. If you do not agree, do not install copy or use the Software. The Software is protected by copyright and other intellectual property laws and is licensed, not sold.    
+
+SCOPE OF RIGHTS:
+You may use, copy, reproduce, and distribute this Software for any non-commercial purpose, subject to the restrictions in this MSR-SSLA. Some purposes which can be non-commercial are teaching, academic research, public demonstrations and personal experimentation. You may also distribute this Software with books or other teaching materials, or publish the Software on websites, that are intended to teach the use of the Software for academic or other non-commercial purposes.
+You may not use or distribute this Software or any derivative works in any form for commercial purposes. Examples of commercial purposes would be running business operations, licensing, leasing, or selling the Software, distributing the Software for use with commercial products, using the Software in the creation or use of commercial products or any other activity which purpose is to procure a commercial gain to you or others.
+If the Software includes source code or data, you may create derivative works of such portions of the Software and distribute the modified Software for non-commercial purposes, as provided herein.
+
+In return, we simply require that you agree: 
+1. That you will not remove any copyright or other notices from the Software.
+2. That if any of the Software is in binary format, you will not attempt to modify such portions of the Software, or to reverse engineer or decompile them, except and only to the extent authorized by applicable law. 
+3. That if you distribute the Software or any derivative works of the Software, you will distribute them under the same terms and conditions as in this license, and you will not grant other rights to the Software or derivative works that are different from those provided by this MSR-SSLA. 
+4. That if you have created derivative works of the Software, and distribute such derivative works, you will cause the modified files to carry prominent notices so that recipients know that they are not receiving the original Software. Such notices must state: (i) that you have changed the Software; and (ii) the date of any changes.
+5. That Microsoft is granted back, without any restrictions or limitations, a non-exclusive, perpetual, irrevocable, royalty-free, assignable and sub-licensable license, to reproduce, publicly perform or display, install, use, modify, distribute, make and have made, sell and transfer your modifications to and/or derivative works of the Software source code or data, for any purpose.  
+6. That any feedback about the Software provided by you to us is voluntarily given, and Microsoft shall be free to use the feedback as it sees fit without obligation or restriction of any kind, even if the feedback is designated by you as confidential. 
+7. THAT THE SOFTWARE COMES "AS IS", WITH NO WARRANTIES. THIS MEANS NO EXPRESS, IMPLIED OR STATUTORY WARRANTY, INCLUDING WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, ANY WARRANTY AGAINST INTERFERENCE WITH YOUR ENJOYMENT OF THE SOFTWARE OR ANY WARRANTY OF TITLE OR NON-INFRINGEMENT. THERE IS NO WARRANTY THAT THIS SOFTWARE WILL FULFILL ANY OF YOUR PARTICULAR PURPOSES OR NEEDS. ALSO, YOU MUST PASS THIS DISCLAIMER ON WHENEVER YOU DISTRIBUTE THE SOFTWARE OR DERIVATIVE WORKS.
+8. THAT NEITHER MICROSOFT NOR ANY CONTRIBUTOR TO THE SOFTWARE WILL BE LIABLE FOR ANY DAMAGES RELATED TO THE SOFTWARE OR THIS MSR-SSLA, INCLUDING DIRECT, INDIRECT, SPECIAL, CONSEQUENTIAL OR INCIDENTAL DAMAGES, TO THE MAXIMUM EXTENT THE LAW PERMITS, NO MATTER WHAT LEGAL THEORY IT IS BASED ON. ALSO, YOU MUST PASS THIS LIMITATION OF LIABILITY ON WHENEVER YOU DISTRIBUTE THE SOFTWARE OR DERIVATIVE WORKS.
+9. That we have no duty of reasonable care or lack of negligence, and we are not obligated to (and will not) provide technical support for the Software.
+10. That if you breach this MSR-SSLA or if you sue anyone over patents that you think may apply to or read on the Software or anyone's use of the Software, this MSR-SSLA (and your license and rights obtained herein) terminate automatically.  Upon any such termination, you shall destroy all of your copies of the Software immediately.  Sections 5, 6, 7, 8, 9, 10, 13 and 14 of this MSR-SSLA shall survive any termination of this MSR-SSLA.
+11. That the patent rights, if any, granted to you in this MSR-SSLA only apply to the Software, not to any derivative works you make.
+12. That the Software may be subject to U.S. export jurisdiction at the time it is licensed to you, and it may be subject to additional export or import laws in other places. You agree to comply with all such laws and regulations that may apply to the Software after delivery of the software to you.
+13. That all rights not expressly granted to you in this MSR-SSLA are reserved.
+14. That this MSR-SSLA shall be construed and controlled by the laws of the State of Washington, USA, without regard to conflicts of law.  If any provision of this MSR-SSLA shall be deemed unenforceable or contrary to law, the rest of this MSR-SSLA shall remain in full effect and interpreted in an enforceable manner that most nearly captures the intent of the original language.
+
+
+Copyright © Microsoft Corporation. All rights reserved.

libs/Math/TRWS/MRFEnergy.h

@@ -0,0 +1,246 @@
+/******************************************************************
+Vladimir Kolmogorov, 2005
+vnk@microsoft.com
+
+(c) Microsoft Corporation. All rights reserved. 
+*******************************************************************/
+
+#ifndef __MRFENERGY_H__
+#define __MRFENERGY_H__
+
+#include "instances.h"
+
+
+// After MRFEnergy is allocated, there are two phases:
+// 1. Energy construction. Only AddNode(), AddNodeData() and AddEdge() may be called.
+// 
+// Any call ZeroMessages(), SetAutomaticOrdering(), Minimize_TRW_S() or Minimize_BP()
+// completes graph construction; MRFEnergy goes to the second phase:
+// 2. Only functions AddNodeData(), ZeroMessages(), Minimize_TRW_S(), Minimize_BP()
+// or GetSolution() may be called. (The last function can be called only after
+// Minimize_TRW_S() or Minimize_BP()).
+
+
+template <class T> class MRFEnergy
+{
+private:
+	struct Node;
+
+public:
+	typedef typename T::Label      Label;
+	typedef typename T::REAL       REAL;
+	typedef typename T::GlobalSize GlobalSize;
+	typedef typename T::LocalSize  LocalSize;
+	typedef typename T::NodeData   NodeData;
+	typedef typename T::EdgeData   EdgeData;
+
+	typedef Node* NodeId;
+	typedef void (*ErrorFunction)(const char* msg);
+
+	// Constructor. Function errorFn is called with an error message, if an error occurs.
+	MRFEnergy(GlobalSize Kglobal, ErrorFunction errorFn = NULL);
+
+	// Destructor.
+	~MRFEnergy();
+
+	//////////////////////////////////////////////////////////
+	//                 Energy construction                  //
+	//////////////////////////////////////////////////////////
+
+	// Adds a node with parameters K and data 
+	// (see the corresponding message*.h file for description).
+	// Note: information in data is copied into internal memory.
+	// Cannot be called after energy construction is completed.
+	NodeId AddNode(LocalSize K, NodeData data);
+
+	// Modifies node parameter for existing node (namely, add information
+	// in data to existing parameter). May be called at any time.
+	// Node i must be NodeId returned by AddNode().
+	void AddNodeData(NodeId i, NodeData data);
+
+	// Adds an edge between i and j. data determines edge parameters
+	// (see the corresponding message*.h file for description).
+	// Note: information in data is copied into internal memory.
+	// Cannot be called after energy construction is completed.
+	void AddEdge(NodeId i, NodeId j, EdgeData data);
+
+	//////////////////////////////////////////////////////////
+	//                Energy construction end               //
+	//////////////////////////////////////////////////////////
+
+	// Clears all messages. Completes energy construction (if not completed yet).
+	void ZeroMessages();
+
+	// Adds to all message entries a value drawn uniformly from [min_value, max_value].
+	// Normally, min_value can be set to 0 (except for TypeBinaryFast, in which case min_value = -max_value)
+	void AddRandomMessages(unsigned int random_seed, REAL min_value, REAL max_value);
+
+	// The algorithm depends on the order of nodes.
+	// By default nodes are processed in the order in which they were added.
+	// The function below permutes this order using certain heuristics.
+	// It may speed up the algorithm if, for example, the original order is random.
+	// 
+	// Completes energy construction.
+	// Cannot be called after energy construction is completed.
+	void SetAutomaticOrdering();
+
+	// The structure below specifies (1) stopping criteria and 
+	// (2) how often to compute solution and print its energy.
+	struct Options
+	{
+		Options()
+		{
+			// default parameters
+			m_eps = -1; // not used
+			m_iterMax = 1000000;
+			m_printIter = 5;     // After 10 iterations start printing the lower bound
+			m_printMinIter = 10; // and the energy every 5 iterations.
+		}
+
+		// stopping criterion
+		REAL		m_eps; // stop if the increase in the lower bound during one iteration is less or equal than m_eps.
+						   // Used only if m_eps >= 0, and only for TRW-S algorithm.
+		int			m_iterMax; // maximum number of iterations
+
+		// Option for printing lower bound and the energy.
+		// Note: computing solution and its energy is slow
+		// (it is comparable to the cost of one iteration).
+		int		m_printIter; // print lower bound and energy every m_printIter iterations
+		int		m_printMinIter; // do not print lower bound and energy before m_printMinIter iterations
+	};
+
+	// Returns number of iterations. Sets lowerBound and energy.
+	// If the user provides array min_marginals, then the code
+	// sets this array accordingly. (The size of the array depends on the type
+	// used. Normally, it's (# nodes)*(# labels). Exception: for TypeBinaryFast it's (# nodes).
+	int Minimize_TRW_S(Options& options, REAL& lowerBound, REAL& energy, REAL* min_marginals = NULL);
+
+	// Returns number of iterations. Sets energy.
+	int Minimize_BP(Options& options, REAL& energy, REAL* min_marginals = NULL);
+
+	// Returns an integer in [0,Ki). Can be called only after Minimize().
+	Label GetSolution(NodeId i);
+
+
+
+
+
+
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////
+	//                   Implementation                     //
+	//////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////
+	//////////////////////////////////////////////////////////
+private:
+
+	typedef typename T::Vector Vector;
+	typedef typename T::Edge   Edge;
+
+	struct MRFEdge;
+	struct MallocBlock;
+
+	ErrorFunction	m_errorFn;
+	MallocBlock*	m_mallocBlockFirst;
+	Node*			m_nodeFirst;
+	Node*			m_nodeLast;
+	int				m_nodeNum;
+	int				m_edgeNum;
+	GlobalSize		m_Kglobal;
+	int				m_vectorMaxSizeInBytes;
+
+	bool			m_isEnergyConstructionCompleted;
+
+	char*			m_buf; // buffer of size m_vectorMaxSizeInBytes 
+					       //              + max(m_vectorMaxSizeInBytes, Edge::GetBufSizeInBytes(m_vectorMaxSizeInBytes))
+
+	void CompleteGraphConstruction(); // nodes and edges cannot be added after calling this function
+	void SetMonotonicTrees();
+
+	REAL ComputeSolutionAndEnergy(); // sets Node::m_solution, returns value of the energy
+
+
+
+	struct Node
+	{
+		int			m_ordering; // unique integer in [0,m_nodeNum-1)
+
+		MRFEdge*	m_firstForward; // first edge going to nodes with greater m_ordering
+		MRFEdge*	m_firstBackward; // first edge going to nodes with smaller m_ordering
+
+		Node*		m_prev; // previous and next
+		Node*		m_next; // nodes according to m_ordering
+
+		Label		m_solution; // integer in [0,m_D.m_K)
+		LocalSize	m_K; // local information about number of labels
+
+		Vector		m_D; // must be the last member in the struct since its size is not fixed
+	};
+
+	struct MRFEdge
+	{
+		MRFEdge*	m_nextForward; // next forward edge with the same tail
+		MRFEdge*	m_nextBackward; // next backward edge with the same head
+		Node*		m_tail;
+		Node*		m_head;
+
+		REAL		m_gammaForward; // = rho_{ij} / rho_{i} where i=m_tail, j=m_head
+		REAL		m_gammaBackward; // = rho_{ij} / rho_{j} where i=m_tail, j=m_head
+
+		Edge		m_message; // must be the last member in the struct since its size is not fixed.
+					           // Stores edge information and either forward or backward message.
+					           // Most of the time it's the backward message; it gets replaced
+					           // by the forward message only temporarily inside Minimize_TRW_S() and Minimize_BP().
+	};
+
+	// Use our own Malloc since 
+	// (a) new in C++ is slow and allocates minimum memory of 64 bytes (in Visual C++)
+	// (b) we want simple (one function) deallocation instead of going through every allocated element
+	struct MallocBlock
+	{
+		static const int minBlockSizeInBytes = 4096 - 3*sizeof(void*);
+		MallocBlock*	m_next;
+		char*			m_current; // first element of available memory in this block
+		char*			m_last; // first element outside of allocated memory for this block
+	};
+	char* Malloc(int bytesNum); 
+};
+
+
+
+
+template <class T> inline char* MRFEnergy<T>::Malloc(int bytesNum)
+{
+	if (!m_mallocBlockFirst || m_mallocBlockFirst->m_current+bytesNum > m_mallocBlockFirst->m_last)
+	{
+		int size = (bytesNum > MallocBlock::minBlockSizeInBytes) ? bytesNum : MallocBlock::minBlockSizeInBytes;
+		MallocBlock* b = (MallocBlock*) new char[sizeof(MallocBlock) + size];
+		if (!b) m_errorFn("Not enough memory");
+		b->m_current = (char*) b + sizeof(MallocBlock);
+		b->m_last    = b->m_current + size;
+
+		b->m_next = m_mallocBlockFirst;
+		m_mallocBlockFirst = b;
+	}
+
+	char* ptr = m_mallocBlockFirst->m_current;
+	m_mallocBlockFirst->m_current += bytesNum;
+	return ptr;
+}
+
+template <class T> inline typename T::Label MRFEnergy<T>::GetSolution(NodeId i)
+{
+	return i->m_solution;
+}
+
+#include "MRFEnergy.inl"
+
+#endif

libs/Math/TRWS/MRFEnergy.inl

@@ -0,0 +1,255 @@
+#include "minimize.inl"
+#include "ordering.inl"
+#include "treeProbabilities.inl"
+
+static void DefaultErrorFn(const char* msg)
+{
+	fprintf(stderr, "%s\n", msg);
+	exit(1);
+}
+
+template <class T> MRFEnergy<T>::MRFEnergy(GlobalSize Kglobal, ErrorFunction errorFn)
+	: m_errorFn(errorFn ? errorFn : DefaultErrorFn),
+	  m_mallocBlockFirst(NULL),
+	  m_nodeFirst(NULL),
+	  m_nodeLast(NULL),
+	  m_nodeNum(0),
+	  m_edgeNum(0),
+	  m_Kglobal(Kglobal),
+	  m_vectorMaxSizeInBytes(0),
+	  m_isEnergyConstructionCompleted(false),
+	  m_buf(NULL)
+{
+}
+
+template <class T> MRFEnergy<T>::~MRFEnergy()
+{
+	while (m_mallocBlockFirst)
+	{
+		MallocBlock* next = m_mallocBlockFirst->m_next;
+		delete m_mallocBlockFirst;
+		m_mallocBlockFirst = next;
+	}
+}
+
+template <class T> typename MRFEnergy<T>::NodeId MRFEnergy<T>::AddNode(LocalSize K, NodeData data)
+{
+	if (m_isEnergyConstructionCompleted)
+	{
+		m_errorFn("Error in AddNode(): graph construction completed - nodes cannot be added");
+	}
+
+	int actualVectorSize = Vector::GetSizeInBytes(m_Kglobal, K);
+	if (actualVectorSize < 0)
+	{
+		m_errorFn("Error in AddNode() (invalid parameter?)");
+	}
+	if (m_vectorMaxSizeInBytes < actualVectorSize)
+	{
+		m_vectorMaxSizeInBytes = actualVectorSize;
+	}
+	int nodeSize = sizeof(Node) - sizeof(Vector) + actualVectorSize;
+	Node* i = (Node *) Malloc(nodeSize);
+
+	i->m_K = K;
+	i->m_D.Initialize(m_Kglobal, K, data);
+
+	i->m_firstForward = NULL;
+	i->m_firstBackward = NULL;
+	i->m_prev = m_nodeLast;
+	if (m_nodeLast)
+	{
+		m_nodeLast->m_next = i;
+	}
+	else
+	{
+		m_nodeFirst = i;
+	}
+	m_nodeLast = i;
+	i->m_next = NULL;
+
+	i->m_ordering = m_nodeNum ++;
+
+	return i;
+}
+
+template <class T> void MRFEnergy<T>::AddNodeData(NodeId i, NodeData data)
+{
+	i->m_D.Add(m_Kglobal, i->m_K, data);
+}
+
+template <class T> void MRFEnergy<T>::AddEdge(NodeId i, NodeId j, EdgeData data)
+{
+	if (m_isEnergyConstructionCompleted)
+	{
+		m_errorFn("Error in AddNode(): graph construction completed - nodes cannot be added");
+	}
+
+	MRFEdge* e;
+
+	int actualEdgeSize = Edge::GetSizeInBytes(m_Kglobal, i->m_K, j->m_K, data);
+	if (actualEdgeSize < 0)
+	{
+		m_errorFn("Error in AddEdge() (invalid parameter?)");
+	}
+	int MRFedgeSize = sizeof(MRFEdge) - sizeof(Edge) + actualEdgeSize;
+	e = (MRFEdge*) Malloc(MRFedgeSize);
+
+	e->m_message.Initialize(m_Kglobal, i->m_K, j->m_K, data, &i->m_D, &j->m_D);
+
+	e->m_tail = i;
+	e->m_nextForward = i->m_firstForward;
+	i->m_firstForward = e;
+
+	e->m_head = j;
+	e->m_nextBackward = j->m_firstBackward;
+	j->m_firstBackward = e;
+
+	m_edgeNum ++;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+
+template <class T> void MRFEnergy<T>::ZeroMessages()
+{
+	Node* i;
+	MRFEdge* e;
+
+	if (!m_isEnergyConstructionCompleted)
+	{
+		CompleteGraphConstruction();
+	}
+
+	for (i=m_nodeFirst; i; i=i->m_next)
+	{
+		for (e=i->m_firstForward; e; e=e->m_nextForward)
+		{
+			e->m_message.GetMessagePtr()->SetZero(m_Kglobal, i->m_K);
+		}
+	}
+}
+
+template <class T> void MRFEnergy<T>::AddRandomMessages(unsigned int random_seed, REAL min_value, REAL max_value)
+{
+	Node* i;
+	MRFEdge* e;
+	int k;
+
+	if (!m_isEnergyConstructionCompleted)
+	{
+		CompleteGraphConstruction();
+	}
+
+	srand(random_seed);
+
+	for (i=m_nodeFirst; i; i=i->m_next)
+	{
+		for (e=i->m_firstForward; e; e=e->m_nextForward)
+		{
+			Vector* M = e->m_message.GetMessagePtr();
+			for (k=0; k<M->GetArraySize(m_Kglobal, i->m_K); k++)
+			{
+				REAL x = (REAL)( min_value + rand()/((double)RAND_MAX) * (max_value - min_value) );
+				x += M->GetArrayValue(m_Kglobal, i->m_K, k);
+				M->SetArrayValue(m_Kglobal, i->m_K, k, x);
+			}
+		}
+	}
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+
+template <class T> void MRFEnergy<T>::CompleteGraphConstruction()
+{
+	Node* i;
+	Node* j;
+	MRFEdge* e;
+	MRFEdge* ePrev;
+
+	if (m_isEnergyConstructionCompleted)
+	{
+		m_errorFn("Fatal error in CompleteGraphConstruction");
+	}
+
+	printf("Completing graph construction... ");
+
+	if (m_buf)
+	{
+		m_errorFn("CompleteGraphConstruction(): fatal error");
+	}
+
+	m_buf = (char *) Malloc(m_vectorMaxSizeInBytes + 
+		( m_vectorMaxSizeInBytes > Edge::GetBufSizeInBytes(m_vectorMaxSizeInBytes) ?
+		  m_vectorMaxSizeInBytes : Edge::GetBufSizeInBytes(m_vectorMaxSizeInBytes) ) );
+
+	// set forward and backward edges properly
+#ifdef _DEBUG
+	int ordering;
+	for (i=m_nodeFirst, ordering=0; i; i=i->m_next, ordering++)
+	{
+		if ( (i->m_ordering != ordering)
+		  || (i->m_ordering == 0 && i->m_prev)
+		  || (i->m_ordering != 0 && i->m_prev->m_ordering != ordering-1) )
+		{
+			m_errorFn("CompleteGraphConstruction(): fatal error (wrong ordering)");
+		}
+	}
+	if (ordering != m_nodeNum)
+	{
+		m_errorFn("CompleteGraphConstruction(): fatal error");
+	}
+#endif
+	for (i=m_nodeFirst; i; i=i->m_next)
+	{
+		i->m_firstBackward = NULL;
+	}
+	for (i=m_nodeFirst; i; i=i->m_next)
+	{
+		ePrev = NULL;
+		for (e=i->m_firstForward; e; )
+		{
+			assert(i == e->m_tail);
+			j = e->m_head;
+
+			if (i->m_ordering < j->m_ordering)
+			{
+				e->m_nextBackward = j->m_firstBackward;
+				j->m_firstBackward = e;
+
+				ePrev = e;
+				e = e->m_nextForward;
+			}
+			else
+			{
+				e->m_message.Swap(m_Kglobal, i->m_K, j->m_K);
+				e->m_tail = j;
+				e->m_head = i;
+
+				MRFEdge* eNext = e->m_nextForward;
+
+				if (ePrev)
+				{
+					ePrev->m_nextForward = e->m_nextForward;
+				}
+				else
+				{
+					i->m_firstForward = e->m_nextForward;
+				}
+
+				e->m_nextForward = j->m_firstForward;
+				j->m_firstForward = e;
+
+				e->m_nextBackward = i->m_firstBackward;
+				i->m_firstBackward = e;
+
+				e = eNext;
+			}
+		}
+	}
+
+	m_isEnergyConstructionCompleted = true;
+
+	// ZeroMessages();
+
+	printf("done\n");
+}

libs/Math/TRWS/README.TXT

@@ -0,0 +1,21 @@
+This software implements two algorithms for minimizing energy functions of the form 
+
+   E(x)   =   \sum_i D_i(x_i)   +   \sum_ij V_ij(x_i,x_j)
+   where x_i are discrete variables.
+
+The two algorithms are max-product belief propagation (BP, Pearl'88) and
+sequential tree-reweighted max-product message passing (TRW-S, Kolmogorov'05).
+
+For example usage look at one of the type*.h files
+(typeBinary.h, typeBinaryFast.h, typePotts.h, typeGeneral.h,
+typeTruncatedLinear.h, typeTruncatedQuadratic.h, typeTruncatedLinear2D.h, typeTruncatedQuadratic2D.h).
+If your energy function does not belong to the classes defined
+in these files but terms V_ij allow fast distance transforms, then
+it should be possible to extend the algorithms to your functions
+using files type*.h as examples.
+
+Written by Vladimir Kolmogorov (vnk@microsoft.com), 2005.
+Tested under Microsoft Visual Studio .NET (Windows)
+and GNU c++ compiler version 2.96 (Red Hat Linux 7.1).
+
+(c) Microsoft Corporation. All rights reserved. 

libs/Math/TRWS/instances.h

@@ -0,0 +1,22 @@
+#ifndef __INSTANCES_H__
+#define __INSTANCES_H__
+
+
+#if defined(_MSC_VER)
+
+// C4661: '...' : no suitable definition provided for explicit template instantiation request
+#pragma warning(disable: 4661)
+
+#endif
+
+#include "typeBinary.h"
+#include "typeBinaryFast.h"
+#include "typePotts.h"
+#include "typeGeneral.h"
+#include "typeTruncatedLinear.h"
+#include "typeTruncatedQuadratic.h"
+#include "typeTruncatedLinear2D.h"
+#include "typeTruncatedQuadratic2D.h"
+
+
+#endif

libs/Math/TRWS/minimize.inl

@@ -0,0 +1,299 @@
+template <class T> int MRFEnergy<T>::Minimize_TRW_S(Options& options, REAL& lowerBound, REAL& energy, REAL* min_marginals)
+{
+	Node* i;
+	Node* j;
+	MRFEdge* e;
+	REAL vMin;
+	int iter;
+	REAL lowerBoundPrev;
+
+	if (!m_isEnergyConstructionCompleted)
+	{
+		CompleteGraphConstruction();
+	}
+
+	printf("TRW_S algorithm\n");
+
+	SetMonotonicTrees();
+
+	Vector* Di = (Vector*) m_buf;
+	void* buf = (void*) (m_buf + m_vectorMaxSizeInBytes);
+
+	iter = 0;
+	bool lastIter = false;
+
+	// main loop
+	for (iter=1; ; iter++)
+	{
+		if (iter >= options.m_iterMax) lastIter = true;
+
+		////////////////////////////////////////////////
+		//                forward pass                //
+		////////////////////////////////////////////////
+		REAL* min_marginals_ptr = min_marginals;
+
+		for (i=m_nodeFirst; i; i=i->m_next)
+		{
+			Di->Copy(m_Kglobal, i->m_K, &i->m_D);
+			for (e=i->m_firstForward; e; e=e->m_nextForward)
+			{
+				Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+			}
+			for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+			{
+				Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+			}
+
+			// normalize Di, update lower bound
+			// vMin = Di->ComputeAndSubtractMin(m_Kglobal, i->m_K); // do not compute lower bound
+			// lowerBound += vMin;                                  // during the forward pass
+
+			// pass messages from i to nodes with higher m_ordering
+			for (e=i->m_firstForward; e; e=e->m_nextForward)
+			{
+				assert(e->m_tail == i);
+				j = e->m_head;
+
+				vMin = e->m_message.UpdateMessage(m_Kglobal, i->m_K, j->m_K, Di, e->m_gammaForward, 0, buf);
+
+				// lowerBound += vMin; // do not compute lower bound during the forward pass
+			}
+
+			if (lastIter && min_marginals)
+			{
+				min_marginals_ptr += Di->GetArraySize(m_Kglobal, i->m_K);
+			}
+		}
+
+		////////////////////////////////////////////////
+		//               backward pass                //
+		////////////////////////////////////////////////
+		lowerBound = 0;
+
+		for (i=m_nodeLast; i; i=i->m_prev)
+		{
+			Di->Copy(m_Kglobal, i->m_K, &i->m_D);
+			for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+			{
+				Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+			}
+			for (e=i->m_firstForward; e; e=e->m_nextForward)
+			{
+				Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+			}
+
+			// normalize Di, update lower bound
+			vMin = Di->ComputeAndSubtractMin(m_Kglobal, i->m_K);
+			lowerBound += vMin;
+
+			// pass messages from i to nodes with smaller m_ordering
+			for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+			{
+				assert(e->m_head == i);
+				j = e->m_tail;
+
+				vMin = e->m_message.UpdateMessage(m_Kglobal, i->m_K, j->m_K, Di, e->m_gammaBackward, 1, buf);
+
+				lowerBound += vMin;
+			}
+
+			if (lastIter && min_marginals)
+			{
+				min_marginals_ptr -= Di->GetArraySize(m_Kglobal, i->m_K);
+				for (int k=0; k<Di->GetArraySize(m_Kglobal, i->m_K); k++) 
+				{
+					min_marginals_ptr[k] = Di->GetArrayValue(m_Kglobal, i->m_K, k);
+				}
+			}
+		}
+
+		////////////////////////////////////////////////
+		//          check stopping criterion          //
+		////////////////////////////////////////////////
+
+		// print lower bound and energy, if necessary
+		if (  lastIter || 
+			( iter>=options.m_printMinIter && 
+			(options.m_printIter<1 || iter%options.m_printIter==0) )
+		)
+		{
+			energy = ComputeSolutionAndEnergy();
+			printf("iter %d: lower bound = %f, energy = %f\n", iter, lowerBound, energy);
+		}
+
+		if (lastIter) break;
+
+		// check convergence of lower bound
+		if (options.m_eps >= 0)
+		{
+			if (iter > 1 && lowerBound - lowerBoundPrev <= options.m_eps)
+			{
+				lastIter = true;
+			}
+			lowerBoundPrev = lowerBound;
+		}
+	}
+
+	return iter;
+}
+
+template <class T> int MRFEnergy<T>::Minimize_BP(Options& options, REAL& energy, REAL* min_marginals)
+{
+	Node* i;
+	Node* j;
+	MRFEdge* e;
+	REAL vMin;
+	int iter;
+
+	if (!m_isEnergyConstructionCompleted)
+	{
+		CompleteGraphConstruction();
+	}
+
+	printf("BP algorithm\n");
+
+	Vector* Di = (Vector*) m_buf;
+	void* buf = (void*) (m_buf + m_vectorMaxSizeInBytes);
+
+	iter = 0;
+	bool lastIter = false;
+
+	// main loop
+	for (iter=1; ; iter++)
+	{
+		if (iter >= options.m_iterMax) lastIter = true;
+
+		////////////////////////////////////////////////
+		//                forward pass                //
+		////////////////////////////////////////////////
+		REAL* min_marginals_ptr = min_marginals;
+
+		for (i=m_nodeFirst; i; i=i->m_next)
+		{
+			Di->Copy(m_Kglobal, i->m_K, &i->m_D);
+			for (e=i->m_firstForward; e; e=e->m_nextForward)
+			{
+				Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+			}
+			for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+			{
+				Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+			}
+
+			// pass messages from i to nodes with higher m_ordering
+			for (e=i->m_firstForward; e; e=e->m_nextForward)
+			{
+				assert(i == e->m_tail);
+				j = e->m_head;
+
+				const REAL gamma = 1;
+
+				e->m_message.UpdateMessage(m_Kglobal, i->m_K, j->m_K, Di, gamma, 0, buf);
+			}
+
+			if (lastIter && min_marginals)
+			{
+				min_marginals_ptr += Di->GetArraySize(m_Kglobal, i->m_K);
+			}
+		}
+
+		////////////////////////////////////////////////
+		//               backward pass                //
+		////////////////////////////////////////////////
+
+		for (i=m_nodeLast; i; i=i->m_prev)
+		{
+			Di->Copy(m_Kglobal, i->m_K, &i->m_D);
+			for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+			{
+				Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+			}
+			for (e=i->m_firstForward; e; e=e->m_nextForward)
+			{
+				Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+			}
+
+			// pass messages from i to nodes with smaller m_ordering
+			for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+			{
+				assert(i == e->m_head);
+				j = e->m_tail;
+
+				const REAL gamma = 1;
+
+				vMin = e->m_message.UpdateMessage(m_Kglobal, i->m_K, j->m_K, Di, gamma, 1, buf);
+			}
+
+			if (lastIter && min_marginals)
+			{
+				min_marginals_ptr -= Di->GetArraySize(m_Kglobal, i->m_K);
+				for (int k=0; k<Di->GetArraySize(m_Kglobal, i->m_K); k++) 
+				{
+					min_marginals_ptr[k] = Di->GetArrayValue(m_Kglobal, i->m_K, k);
+				}
+			}
+		}
+
+		////////////////////////////////////////////////
+		//          check stopping criterion          //
+		////////////////////////////////////////////////
+
+		// print energy, if necessary
+		if ( lastIter || 
+			( iter>=options.m_printMinIter && 
+			(options.m_printIter<1 || iter%options.m_printIter==0) )
+		)
+		{
+			energy = ComputeSolutionAndEnergy();
+			printf("iter %d: energy = %f\n", iter, energy);
+		}
+
+		// if finishFlag==true terminate
+		if (lastIter) break;
+	}
+
+	return iter;
+}
+
+template <class T> typename T::REAL MRFEnergy<T>::ComputeSolutionAndEnergy()
+{
+	Node* i;
+	Node* j;
+	MRFEdge* e;
+	REAL E = 0;
+
+	Vector* DiBackward = (Vector*) m_buf; // cost of backward edges plus Di at the node
+	Vector* Di = (Vector*) (m_buf + m_vectorMaxSizeInBytes); // all edges plus Di at the node
+
+	for (i=m_nodeFirst; i; i=i->m_next)
+	{
+		// Set Ebackward[ki] to be the sum of V(ki,j->m_solution) for backward edges (i,j).
+		// Set Di[ki] to be the value of the energy corresponding to
+		// part of the graph considered so far, assuming that nodes u
+		// in this subgraph are fixed to u->m_solution
+
+		DiBackward->Copy(m_Kglobal, i->m_K, &i->m_D);
+		for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+		{
+			assert(i == e->m_head);
+			j = e->m_tail;
+
+			e->m_message.AddColumn(m_Kglobal, j->m_K, i->m_K, j->m_solution, DiBackward, 0);
+		}
+
+		// add forward edges
+		Di->Copy(m_Kglobal, i->m_K, DiBackward);
+
+		for (e=i->m_firstForward; e; e=e->m_nextForward)
+		{
+			Di->Add(m_Kglobal, i->m_K, e->m_message.GetMessagePtr());
+		}
+
+		Di->ComputeMin(m_Kglobal, i->m_K, i->m_solution);
+
+		// update energy
+		E += DiBackward->GetValue(m_Kglobal, i->m_K, i->m_solution);
+	}
+
+	return E;
+}

libs/Math/TRWS/ordering.inl

@@ -0,0 +1,151 @@
+template <class T> void MRFEnergy<T>::SetAutomaticOrdering()
+{
+	int dMin;
+	Node* i;
+	Node* iMin;
+	Node* list;
+	Node* listBoundary;
+	MRFEdge* e;
+
+	if (m_isEnergyConstructionCompleted)
+	{
+		m_errorFn("Error in SetAutomaticOrdering(): function cannot be called after graph construction is completed");
+	}
+
+	printf("Setting automatic ordering... ");
+
+	list = m_nodeFirst;
+	listBoundary = NULL;
+	m_nodeFirst = m_nodeLast = NULL;
+	for (i=list; i; i=i->m_next)
+	{
+		i->m_ordering = 2*m_nodeNum; // will contain remaining degree mod m_nodeNum (i.e. number of edges connecting to nodes in 'listBoundary' and 'list')
+		                             // if i->m_ordering \in [2*m_nodeNum;  3*m_nodeNum) - not assigned yet, belongs to 'list'
+		                             // if i->m_ordering \in [m_nodeNum;    2*m_nodeNum) - not assigned yet, belongs to 'listBoundary'
+		                             // if i->m_ordering \in [0;            m_nodeNum  ) - assigned, belongs to 'm_nodeFirst'
+		for (e=i->m_firstForward; e; e=e->m_nextForward)
+		{
+			i->m_ordering ++;
+		}
+		for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+		{
+			i->m_ordering ++;
+		}
+	}
+
+	while (list)
+	{
+		// find node with the smallest remaining degree in list
+		dMin = m_nodeNum;
+		for (i=list; i; i=i->m_next)
+		{
+			assert(i->m_ordering >= 2*m_nodeNum);
+			if (dMin > i->m_ordering - 2*m_nodeNum)
+			{
+				dMin = i->m_ordering - 2*m_nodeNum;
+				iMin = i;
+			}
+		}
+		i = iMin;
+
+		// remove i from list
+		if (i->m_prev) i->m_prev->m_next = i->m_next;
+		else           list              = i->m_next;
+		if (i->m_next) i->m_next->m_prev = i->m_prev;
+
+		// add i to listBoundary
+		listBoundary = i;
+		i->m_prev = NULL;
+		i->m_next = NULL;
+		i->m_ordering -= m_nodeNum;
+
+		while (listBoundary)
+		{
+			// find node with the smallest remaining degree in listBoundary
+			dMin = m_nodeNum;
+			for (i=listBoundary; i; i=i->m_next)
+			{
+				assert(i->m_ordering >= m_nodeNum && i->m_ordering < 2*m_nodeNum);
+				if (dMin > i->m_ordering - m_nodeNum)
+				{
+					dMin = i->m_ordering - m_nodeNum;
+					iMin = i;
+				}
+			}
+			i = iMin;
+
+			// remove i from listBoundary
+			if (i->m_prev) i->m_prev->m_next = i->m_next;
+			else           listBoundary      = i->m_next;
+			if (i->m_next) i->m_next->m_prev = i->m_prev;
+
+			// add i to m_nodeFirst
+			if (m_nodeLast)
+			{
+				m_nodeLast->m_next = i;
+				i->m_ordering = m_nodeLast->m_ordering + 1;
+			}
+			else
+			{
+				m_nodeFirst = i;
+				i->m_ordering = 0;
+			}
+			i->m_prev = m_nodeLast;
+			m_nodeLast = i;
+			i->m_next = NULL;
+
+			// process neighbors of i=m_nodeLast: decrease their remaining degree,
+			// put them into listBoundary (if they are in list)
+			for (e=m_nodeLast->m_firstForward; e; e=e->m_nextForward)
+			{
+				assert(m_nodeLast == e->m_tail);
+				i = e->m_head;
+				if (i->m_ordering >= m_nodeNum)
+				{
+					i->m_ordering --; // decrease remaining degree of i
+					if (i->m_ordering >= 2*m_nodeNum)
+					{
+						// remove i from list
+						if (i->m_prev) i->m_prev->m_next = i->m_next;
+						else           list              = i->m_next;
+						if (i->m_next) i->m_next->m_prev = i->m_prev;
+
+						// add i to listBoundary
+						if (listBoundary) listBoundary->m_prev = i;
+						i->m_prev = NULL;
+						i->m_next = listBoundary;
+						listBoundary = i;
+						i->m_ordering -= m_nodeNum;
+					}
+				}
+			}
+			for (e=m_nodeLast->m_firstBackward; e; e=e->m_nextBackward)
+			{
+				assert(m_nodeLast == e->m_head);
+				i = e->m_tail;
+				if (i->m_ordering >= m_nodeNum)
+				{
+					i->m_ordering --; // decrease remaining degree of i
+					if (i->m_ordering >= 2*m_nodeNum)
+					{
+						// remove i from list
+						if (i->m_prev) i->m_prev->m_next = i->m_next;
+						else           list              = i->m_next;
+						if (i->m_next) i->m_next->m_prev = i->m_prev;
+
+						// add i to listBoundary
+						if (listBoundary) listBoundary->m_prev = i;
+						i->m_prev = NULL;
+						i->m_next = listBoundary;
+						listBoundary = i;
+						i->m_ordering -= m_nodeNum;
+					}
+				}
+			}
+		}
+	}
+
+	printf("done\n");
+
+	CompleteGraphConstruction();
+}

libs/Math/TRWS/treeProbabilities.inl

@@ -0,0 +1,36 @@
+template <class T> void MRFEnergy<T>::SetMonotonicTrees()
+{
+	Node* i;
+	MRFEdge* e;
+
+	if (!m_isEnergyConstructionCompleted)
+	{
+		CompleteGraphConstruction();
+	}
+
+	for (i=m_nodeFirst; i; i=i->m_next)
+	{
+		REAL mu;
+
+		int nForward = 0, nBackward = 0;
+		for (e=i->m_firstForward; e; e=e->m_nextForward)
+		{
+			nForward ++;
+		}
+		for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+		{
+			nBackward ++;
+		}
+		int ni = (nForward > nBackward) ? nForward : nBackward;
+
+		mu = (REAL)1 / ni;
+		for (e=i->m_firstBackward; e; e=e->m_nextBackward)
+		{
+			e->m_gammaBackward = mu;
+		}
+		for (e=i->m_firstForward; e; e=e->m_nextForward)
+		{
+			e->m_gammaForward = mu;
+		}
+	}
+}

libs/Math/TRWS/typeBinary.h

@@ -0,0 +1,392 @@
+/******************************************************************
+typeBinary.h
+
+Energy function with binary labels:
+   E(x)   =   \sum_i D_i(x_i)   +   \sum_ij V_ij(x_i,x_j)
+   where x_i \in {0,1}.
+
+Example usage:
+
+Minimize function E(x,y) = Dx(x) + Dy(y) + V(x,y) where 
+  x,y \in {0,1},
+  Dx(0) = 0, Dx(1) = 1,
+  Dy(0) = 2, Dy(1) = 3,
+  V(0,0) = 4, V(0,1) = 5,
+  V(1,0) = 6, V(1,1) = 7,
+  [.] is 1 if it's argument is true, and 0 otherwise.
+
+
+
+#include <stdio.h>
+#include "MRFEnergy.h"
+
+void testBinary()
+{
+	MRFEnergy<TypeBinary>* mrf;
+	MRFEnergy<TypeBinary>::NodeId* nodes;
+	MRFEnergy<TypeBinary>::Options options;
+	TypeBinary::REAL energy, lowerBound;
+
+	const int nodeNum = 2; // number of nodes
+	int x, y;
+
+	mrf = new MRFEnergy<TypeBinary>(TypeBinary::GlobalSize(K));
+	nodes = new MRFEnergy<TypeBinary>::NodeId[nodeNum];
+
+	// construct energy
+	nodes[0] = mrf->AddNode(TypeBinary::LocalSize(), TypeBinary::NodeData(0, 1));
+	nodes[1] = mrf->AddNode(TypeBinary::LocalSize(), TypeBinary::NodeData(2, 3));
+	mrf->AddEdge(nodes[0], nodes[1], TypeBinary::EdgeData(4, 5, 6, 7));
+
+	// 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\n", x, y);
+
+	// done
+	delete nodes;
+	delete mrf;
+}
+
+*******************************************************************/
+
+
+
+
+
+
+#ifndef __TYPEBINARY_H__
+#define __TYPEBINARY_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypeBinary
+{
+private:
+	struct Vector; // node parameters and messages
+	struct Edge; // stores edge information and either forward or backward message
+
+public:
+	// types declarations
+	typedef int Label;
+	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 // always 2 labels - no need to store it
+	{
+	};
+
+	struct LocalSize // always 2 labels - no need to store it
+	{
+	};
+
+	struct NodeData
+	{
+		NodeData(REAL D0, REAL D1);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_data[2];
+	};
+
+	struct EdgeData
+	{
+		EdgeData(REAL V00, REAL V01, REAL V10, REAL V11);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_V[2][2];
+	};
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	////////////////////////// Visible only to MRFEnergy /////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+
+private:
+friend class MRFEnergy<TypeBinary>;
+
+	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 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[2];
+	};
+
+	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:
+		// edge information
+		REAL		m_lambdaIsing;
+
+		// message
+		Vector		m_message;
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypeBinary::NodeData::NodeData(REAL D0, REAL D1)
+{
+	m_data[0] = D0;
+	m_data[1] = D1;
+}
+
+inline TypeBinary::EdgeData::EdgeData(REAL V00, REAL V01, REAL V10, REAL V11)
+{
+	m_V[0][0] = V00;
+	m_V[0][1] = V01;
+	m_V[1][0] = V10;
+	m_V[1][1] = V11;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypeBinary::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	return sizeof(Vector);
+}
+inline void TypeBinary::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	m_data[0] = data.m_data[0];
+	m_data[1] = data.m_data[1];
+}
+
+inline void TypeBinary::Vector::Add(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	m_data[0] += data.m_data[0];
+	m_data[1] += data.m_data[1];
+}
+
+inline void TypeBinary::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	m_data[0] = 0;
+	m_data[1] = 0;
+}
+
+inline void TypeBinary::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	m_data[0] = V->m_data[0];
+	m_data[1] = V->m_data[1];
+}
+
+inline void TypeBinary::Vector::Add(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	m_data[0] += V->m_data[0];
+	m_data[1] += V->m_data[1];
+}
+
+inline TypeBinary::REAL TypeBinary::Vector::GetValue(GlobalSize Kglobal, LocalSize K, Label k)
+{
+	assert(k>=0 && k<2);
+	return m_data[k];
+}
+
+inline TypeBinary::REAL TypeBinary::Vector::ComputeMin(GlobalSize Kglobal, LocalSize K, Label& kMin)
+{
+	kMin = (m_data[0] <= m_data[1]) ? 0 : 1;
+	return m_data[kMin];
+}
+
+inline TypeBinary::REAL TypeBinary::Vector::ComputeAndSubtractMin(GlobalSize Kglobal, LocalSize K)
+{
+	REAL vMin;
+
+	if (m_data[0] <= m_data[1])
+	{
+		vMin = m_data[0];
+		m_data[0] = 0;
+		m_data[1] -= vMin;
+	}
+	else
+	{
+		vMin = m_data[1];
+		m_data[1] = 0;
+		m_data[0] -= vMin;
+	}
+
+	return vMin;
+}
+
+inline int TypeBinary::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return 2;
+}
+
+inline TypeBinary::REAL TypeBinary::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k>=0 && k<2);
+	return m_data[k];
+}
+
+inline void TypeBinary::Vector::SetArrayValue(GlobalSize Kglobal, LocalSize K, int k, REAL x)
+{
+	assert(k>=0 && k<2);
+	m_data[k] = x;
+}
+
+///////////////////// EdgeDataAndMessage implementation /////////////////////////
+
+inline int TypeBinary::Edge::GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data)
+{
+	return sizeof(Edge);
+}
+
+inline int TypeBinary::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	return 0;
+}
+
+inline void TypeBinary::Edge::Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj)
+{
+	// V00 V01  = A B  = A A  +  0.5 * ( 0   0   + 0 P-Q  +  0   P+Q  )
+	// V10 V11    C D    D D             Q-P Q-P   0 P-Q     P+Q 0
+    // where P=B-A, Q=C-D
+	Di->m_data[0] += data.m_V[0][0];
+	Di->m_data[1] += data.m_V[1][1];
+	REAL P = data.m_V[0][1] - data.m_V[0][0], Q = data.m_V[1][0] - data.m_V[1][1];
+	REAL halfPplusQ = (REAL)0.5 * (P + Q);
+	REAL halfPminusQ = (REAL)0.5 * (P - Q);
+	Di->m_data[1] += -halfPminusQ;
+	Dj->m_data[1] += halfPminusQ;
+	m_lambdaIsing = halfPplusQ;
+	m_message.m_data[0] = 0;
+	m_message.m_data[1] = halfPminusQ + (data.m_V[0][0] - data.m_V[1][1]);
+}
+
+inline TypeBinary::Vector* TypeBinary::Edge::GetMessagePtr()
+{
+	return &m_message;
+}
+
+inline void TypeBinary::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+}
+
+inline TypeBinary::REAL TypeBinary::Edge::UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* buf)
+{
+	REAL data[2], vMin;
+
+	data[0] = gamma*source->m_data[0] - m_message.m_data[0];
+	data[1] = gamma*source->m_data[1] - m_message.m_data[1];
+	m_message.m_data[0] = (data[0] < data[1] + m_lambdaIsing) ? data[0] : data[1] + m_lambdaIsing;
+	m_message.m_data[1] = (data[1] < data[0] + m_lambdaIsing) ? data[1] : data[0] + m_lambdaIsing;
+
+	vMin = (m_message.m_data[0] < m_message.m_data[1]) ? m_message.m_data[0] : m_message.m_data[1];
+	m_message.m_data[0] -= vMin;
+	m_message.m_data[1] -= vMin;
+	return vMin;
+}
+
+
+
+
+inline void TypeBinary::Edge::AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir)
+{
+	dest->m_data[1-ksource] += m_lambdaIsing;
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+#endif

libs/Math/TRWS/typeBinaryFast.h

@@ -0,0 +1,309 @@
+/******************************************************************
+typeBinaryFast.h
+
+Same energy and interface as in typeBinary.h.
+Faster implementation; however, values of lower bound and energy are computed
+incorrectly. (Note that the solution is the same as with typeBinary.h).
+
+For example usage, see typeBinary.h (only replace TypeBinary with TypeBinaryFast).
+*******************************************************************/
+
+
+
+
+
+
+#ifndef __TYPEBINARYFAST_H__
+#define __TYPEBINARYFAST_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypeBinaryFast
+{
+private:
+	struct Vector; // node parameters and messages
+	struct Edge; // stores edge information and either forward or backward message
+
+public:
+	// types declarations
+	typedef int Label;
+	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 // always 2 labels - no need to store it
+	{
+	};
+
+	struct LocalSize // always 2 labels - no need to store it
+	{
+	};
+
+	struct NodeData
+	{
+		NodeData(REAL D0, REAL D1);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_data[2];
+	};
+
+	struct EdgeData
+	{
+		EdgeData(REAL V00, REAL V01, REAL V10, REAL V11);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_V[2][2];
+	};
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	////////////////////////// Visible only to MRFEnergy /////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+
+private:
+friend class MRFEnergy<TypeBinaryFast>;
+
+	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 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; // = D[1] - D[0]
+	};
+
+	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:
+		// edge information
+		REAL		m_lambdaIsing;
+
+		// message
+		Vector		m_message;
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypeBinaryFast::NodeData::NodeData(REAL D0, REAL D1)
+{
+	m_data[0] = D0;
+	m_data[1] = D1;
+}
+
+inline TypeBinaryFast::EdgeData::EdgeData(REAL V00, REAL V01, REAL V10, REAL V11)
+{
+	m_V[0][0] = V00;
+	m_V[0][1] = V01;
+	m_V[1][0] = V10;
+	m_V[1][1] = V11;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypeBinaryFast::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	return sizeof(Vector);
+}
+inline void TypeBinaryFast::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	m_data = data.m_data[1] - data.m_data[0];
+}
+
+inline void TypeBinaryFast::Vector::Add(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	m_data += data.m_data[1] - data.m_data[0];
+}
+
+inline void TypeBinaryFast::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	m_data = 0;
+}
+
+inline void TypeBinaryFast::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	m_data = V->m_data;
+}
+
+inline void TypeBinaryFast::Vector::Add(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	m_data += V->m_data;
+}
+
+inline TypeBinaryFast::REAL TypeBinaryFast::Vector::GetValue(GlobalSize Kglobal, LocalSize K, Label k)
+{
+	assert(k>=0 && k<2);
+	return (k == 0) ? 0 : m_data;
+}
+
+inline TypeBinaryFast::REAL TypeBinaryFast::Vector::ComputeMin(GlobalSize Kglobal, LocalSize K, Label& kMin)
+{
+	kMin = (m_data >= 0) ? 0 : 1;
+	return 0;
+}
+
+inline TypeBinaryFast::REAL TypeBinaryFast::Vector::ComputeAndSubtractMin(GlobalSize Kglobal, LocalSize K)
+{
+	return 0;
+}
+
+inline int TypeBinaryFast::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return 1;
+}
+
+inline TypeBinaryFast::REAL TypeBinaryFast::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k==0);
+	return m_data;
+}
+
+inline void TypeBinaryFast::Vector::SetArrayValue(GlobalSize Kglobal, LocalSize K, int k, REAL x)
+{
+	assert(k==0);
+	m_data = x;
+}
+
+///////////////////// EdgeDataAndMessage implementation /////////////////////////
+
+inline int TypeBinaryFast::Edge::GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data)
+{
+	return sizeof(Edge);
+}
+
+inline int TypeBinaryFast::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	return 0;
+}
+
+inline void TypeBinaryFast::Edge::Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj)
+{
+	// V00 V01  = A B  = A A  +  0.5 * ( 0   0   + 0 P-Q  +  0   P+Q  )
+	// V10 V11    C D    D D             Q-P Q-P   0 P-Q     P+Q 0
+    // where P=B-A, Q=C-D
+	Di->m_data += data.m_V[1][1] - data.m_V[0][0];
+	REAL P = data.m_V[0][1] - data.m_V[0][0], Q = data.m_V[1][0] - data.m_V[1][1];
+	REAL halfPplusQ = (REAL)0.5 * (P + Q);
+	REAL halfPminusQ = (REAL)0.5 * (P - Q);
+	Di->m_data += -halfPminusQ;
+	Dj->m_data += halfPminusQ;
+	m_lambdaIsing = halfPplusQ;
+	m_message.m_data = halfPminusQ + (data.m_V[0][0] - data.m_V[1][1]);
+}
+
+inline TypeBinaryFast::Vector* TypeBinaryFast::Edge::GetMessagePtr()
+{
+	return &m_message;
+}
+
+inline void TypeBinaryFast::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+}
+
+inline TypeBinaryFast::REAL TypeBinaryFast::Edge::UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* buf)
+{
+	REAL s = gamma*source->m_data - m_message.m_data;
+
+	if (m_lambdaIsing >= 0)
+	{
+		if      (s >= m_lambdaIsing)  m_message.m_data = m_lambdaIsing;
+		else if (s <= -m_lambdaIsing) m_message.m_data = -m_lambdaIsing;
+		else                          m_message.m_data = s;
+	}
+	else
+	{
+		if      (s >= -m_lambdaIsing) m_message.m_data = m_lambdaIsing;
+		else if (s <= m_lambdaIsing)  m_message.m_data = -m_lambdaIsing;
+		else                          m_message.m_data = -s;
+	}
+
+	return 0;
+}
+
+inline void TypeBinaryFast::Edge::AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir)
+{
+	dest->m_data += (ksource == 0) ? m_lambdaIsing : -m_lambdaIsing;
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+#endif

libs/Math/TRWS/typeGeneral.h

@@ -0,0 +1,614 @@
+/******************************************************************
+typeGeneral.h
+
+Energy function with general interactions:
+   E(x)   =   \sum_i D_i(x_i)   +   \sum_ij V_ij(x_i,x_j)
+   where x_i \in {0, 1, ..., Ki-1}
+   V_ij(ki, kj) are given either as matrices Ki*Kj, or as Potts terms 
+   (V_ij(ki, kj) = 0 if ki==kj, and lambda_ij otherwise, with non-negative lambda_ij).
+
+   Inefficient! If possible, use other type*.h files.
+
+
+Example usage:
+
+Minimize function E(x,y) = Dx(x) + Dy(y) + lambda*[x != y] + Dz(z) + V(y,z) where 
+  x,y \in {0,1,2}, z \in {0,1}
+  Dx(0) = 0, Dx(1) = 1, Dx(2) = 2,
+  Dy(0) = 3, Dy(1) = 4, Dy(2) = 5,
+  lambda = 6,
+  [.] is 1 if it's argument is true, and 0 otherwise.
+  Dz(0) = 7, Dz(1) = 8,
+  V(y,z) = y*y + z
+
+
+
+#include <stdio.h>
+#include "MRFEnergy.h"
+
+void testGeneral()
+{
+	MRFEnergy<TypeGeneral>* mrf;
+	MRFEnergy<TypeGeneral>::NodeId* nodes;
+	MRFEnergy<TypeGeneral>::Options options;
+	TypeGeneral::REAL energy, lowerBound;
+
+	const int nodeNum = 3; // number of nodes
+	TypeGeneral::REAL Dx[3];
+	TypeGeneral::REAL Dy[3];
+	TypeGeneral::REAL Dz[2];
+	TypeGeneral::REAL V[3*2];
+	int x, y, z;
+
+	mrf = new MRFEnergy<TypeGeneral>(TypeGeneral::GlobalSize());
+	nodes = new MRFEnergy<TypeGeneral>::NodeId[nodeNum];
+
+	// construct energy
+	Dx[0] = 0; Dx[1] = 1; Dx[2] = 2;
+	nodes[0] = mrf->AddNode(TypeGeneral::LocalSize(3), TypeGeneral::NodeData(Dx));
+	Dy[0] = 3; Dy[1] = 4; Dy[2] = 5;
+	nodes[1] = mrf->AddNode(TypeGeneral::LocalSize(3), TypeGeneral::NodeData(Dy));
+	mrf->AddEdge(nodes[0], nodes[1], TypeGeneral::EdgeData(TypeGeneral::POTTS, 6));
+	Dz[0] = 7; Dz[1] = 8;
+	nodes[2] = mrf->AddNode(TypeGeneral::LocalSize(2), TypeGeneral::NodeData(Dz));
+	for (y=0; y<3; y++)
+	{
+		for (z=0; z<2; z++)
+		{
+			V[y + z*3] = y*y + z;
+		}
+	}
+	mrf->AddEdge(nodes[1], nodes[2], TypeGeneral::EdgeData(TypeGeneral::GENERAL, V));
+
+	// 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\n", x, y);
+
+	// done
+	delete nodes;
+	delete mrf;
+}
+
+*******************************************************************/
+
+
+
+
+
+
+
+
+
+
+
+
+#ifndef __TYPEGENERAL_H__
+#define __TYPEGENERAL_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypeGeneral
+{
+private:
+	struct Vector; // node parameters and messages
+	struct Edge; // stores edge information and either forward or backward message
+
+public:
+	typedef enum
+	{
+		GENERAL, // edge information is stored as Ki*Kj matrix. Inefficient!
+		POTTS    // edge information is stored as one number (lambdaPotts).
+	} Type;
+
+	// types declarations
+	typedef int Label;
+	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
+	{
+	};
+
+	struct LocalSize // number of labels is stored at MRFEnergy::m_Kglobal
+	{
+		LocalSize(int K);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		int		m_K; // number of labels
+	};
+
+	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(Type type, REAL lambdaPotts); // type must be POTTS
+		EdgeData(Type type, REAL* data); // type must be GENERAL. data = pointer to array of size Ki*Kj
+		                                 // such that V(ki,kj) = data[ki + Ki*kj]
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		Type		m_type;
+		union
+		{
+			REAL	m_lambdaPotts;
+			REAL*	m_dataGeneral;
+		};
+	};
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	////////////////////////// Visible only to MRFEnergy /////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+
+private:
+friend class MRFEnergy<TypeGeneral>;
+
+	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 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);
+
+	protected:
+
+		Type		m_type;
+
+		// message
+		Vector*		m_message;
+	};
+
+	struct EdgePotts : Edge
+	{
+	private:
+	friend struct Edge;
+		REAL	m_lambdaPotts;
+	};
+
+	struct EdgeGeneral : Edge
+	{
+	private:
+	friend struct Edge;
+		int		m_dir; // 0 if Swap() was called even number of times, 1 otherwise
+		REAL	m_data[1]; // array of size Ki*Kj
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+inline TypeGeneral::LocalSize::LocalSize(int K)
+{
+	m_K = K;
+}
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypeGeneral::NodeData::NodeData(REAL* data)
+{
+	m_data = data;
+}
+
+inline TypeGeneral::EdgeData::EdgeData(Type type, REAL lambdaPotts)
+{
+	assert(type == POTTS);
+	m_type = type;
+	m_lambdaPotts = lambdaPotts;
+}
+
+inline TypeGeneral::EdgeData::EdgeData(Type type, REAL* data)
+{
+	assert(type == GENERAL);
+	m_type = type;
+	m_dataGeneral = data;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypeGeneral::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	if (K.m_K < 1)
+	{
+		return -1;
+	}
+	return K.m_K*sizeof(REAL);
+}
+inline void TypeGeneral::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	memcpy(m_data, data.m_data, K.m_K*sizeof(REAL));
+}
+
+inline void TypeGeneral::Vector::Add(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	for (int k=0; k<K.m_K; k++)
+	{
+		m_data[k] += data.m_data[k];
+	}
+}
+
+inline void TypeGeneral::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	memset(m_data, 0, K.m_K*sizeof(REAL));
+}
+
+inline void TypeGeneral::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	memcpy(m_data, V->m_data, K.m_K*sizeof(REAL));
+}
+
+inline void TypeGeneral::Vector::Add(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	for (int k=0; k<K.m_K; k++)
+	{
+		m_data[k] += V->m_data[k];
+	}
+}
+
+inline TypeGeneral::REAL TypeGeneral::Vector::GetValue(GlobalSize Kglobal, LocalSize K, Label k)
+{
+	assert(k>=0 && k<K.m_K);
+	return m_data[k];
+}
+
+inline TypeGeneral::REAL TypeGeneral::Vector::ComputeMin(GlobalSize Kglobal, LocalSize K, Label& kMin)
+{
+	REAL vMin = m_data[0];
+	kMin = 0;
+	for (int k=1; k<K.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+			kMin = k;
+		}
+	}
+
+	return vMin;
+}
+
+inline TypeGeneral::REAL TypeGeneral::Vector::ComputeAndSubtractMin(GlobalSize Kglobal, LocalSize K)
+{
+	REAL vMin = m_data[0];
+	for (int k=1; k<K.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+		}
+	}
+	for (int k=0; k<K.m_K; k++)
+	{
+		m_data[k] -= vMin;
+	}
+
+	return vMin;
+}
+
+inline int TypeGeneral::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return K.m_K;
+}
+
+inline TypeGeneral::REAL TypeGeneral::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k>=0 && k<K.m_K);
+	return m_data[k];
+}
+
+inline void TypeGeneral::Vector::SetArrayValue(GlobalSize Kglobal, LocalSize K, int k, REAL x)
+{
+	assert(k>=0 && k<K.m_K);
+	m_data[k] = x;
+}
+
+///////////////////// EdgeDataAndMessage implementation /////////////////////////
+
+inline int TypeGeneral::Edge::GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data)
+{
+	int messageSizeInBytes = ((Ki.m_K > Kj.m_K) ? Ki.m_K : Kj.m_K)*sizeof(REAL);
+
+	switch (data.m_type)
+	{
+		case POTTS:
+			if (Ki.m_K != Kj.m_K || data.m_lambdaPotts < 0)
+			{
+				return -1;
+			}
+			return sizeof(EdgePotts) + messageSizeInBytes;
+		case GENERAL:
+			return sizeof(EdgeGeneral) - sizeof(REAL) + Ki.m_K*Kj.m_K*sizeof(REAL) + messageSizeInBytes;
+		default:
+			return -1;
+	}
+}
+
+inline int TypeGeneral::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	return vectorMaxSizeInBytes;
+}
+
+inline void TypeGeneral::Edge::Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj)
+{
+	m_type = data.m_type;
+
+	switch (m_type)
+	{
+		case POTTS:
+			((EdgePotts*)this)->m_lambdaPotts = data.m_lambdaPotts;
+			m_message = (Vector*)((char*)this + sizeof(EdgePotts));
+			break;
+		case GENERAL:
+			((EdgeGeneral*)this)->m_dir = 0;
+			memcpy(((EdgeGeneral*)this)->m_data, data.m_dataGeneral, Ki.m_K*Kj.m_K*sizeof(REAL));
+			m_message = (Vector*)((char*)this + sizeof(EdgeGeneral) - sizeof(REAL) + Ki.m_K*Kj.m_K*sizeof(REAL));
+			break;
+		default:
+			assert(0);
+	}
+
+	memset(m_message->m_data, 0, ((Ki.m_K > Kj.m_K) ? Ki.m_K : Kj.m_K)*sizeof(REAL));
+}
+
+inline TypeGeneral::Vector* TypeGeneral::Edge::GetMessagePtr()
+{
+	return m_message;
+}
+
+inline void TypeGeneral::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+	if (m_type == GENERAL)
+	{
+		((EdgeGeneral*)this)->m_dir = 1 - ((EdgeGeneral*)this)->m_dir;
+	}
+}
+
+inline TypeGeneral::REAL TypeGeneral::Edge::UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* _buf)
+{
+	Vector* buf = (Vector*) _buf;
+	REAL vMin;
+
+	if (m_type == POTTS)
+	{
+		assert(Ksource.m_K == Kdest.m_K);
+
+		int k;
+
+		m_message->m_data[0] = gamma*source->m_data[0] - m_message->m_data[0];
+		vMin = m_message->m_data[0];
+
+		for (k=1; k<Ksource.m_K; k++)
+		{
+			m_message->m_data[k] = gamma*source->m_data[k] - m_message->m_data[k];
+			vMin = buf->m_data[0];
+			if (vMin > m_message->m_data[k])
+			{
+				vMin = m_message->m_data[k];
+			}
+		}
+
+		for (k=0; k<Ksource.m_K; k++)
+		{
+			m_message->m_data[k] -= vMin;
+			if (m_message->m_data[k] > ((EdgePotts*)this)->m_lambdaPotts)
+			{
+				m_message->m_data[k] = ((EdgePotts*)this)->m_lambdaPotts;
+			}
+		}
+	}
+	else if (m_type == GENERAL)
+	{
+		int ksource, kdest;
+		REAL* data = ((EdgeGeneral*)this)->m_data;
+
+		for (ksource=0; ksource<Ksource.m_K; ksource++)
+		{
+			buf->m_data[ksource] = gamma*source->m_data[ksource] - m_message->m_data[ksource];
+		}
+
+		if (dir == ((EdgeGeneral*)this)->m_dir)
+		{
+			for (kdest=0; kdest<Kdest.m_K; kdest++)
+			{
+				vMin = buf->m_data[0] + data[0 + kdest*Ksource.m_K];
+				for (ksource=1; ksource<Ksource.m_K; ksource++)
+				{
+					if (vMin > buf->m_data[ksource] + data[ksource + kdest*Ksource.m_K])
+					{
+						vMin = buf->m_data[ksource] + data[ksource + kdest*Ksource.m_K];
+					}
+				}
+				m_message->m_data[kdest] = vMin;
+			}
+		}
+		else
+		{
+			for (kdest=0; kdest<Kdest.m_K; kdest++)
+			{
+				vMin = buf->m_data[0] + data[kdest + 0*Kdest.m_K];
+				for (ksource=1; ksource<Ksource.m_K; ksource++)
+				{
+					if (vMin > buf->m_data[ksource] + data[kdest + ksource*Kdest.m_K])
+					{
+						vMin = buf->m_data[ksource] + data[kdest + ksource*Kdest.m_K];
+					}
+				}
+				m_message->m_data[kdest] = vMin;
+			}
+		}
+
+		vMin = m_message->m_data[0];
+		for (kdest=1; kdest<Kdest.m_K; kdest++)
+		{
+			if (vMin > m_message->m_data[kdest])
+			{
+				vMin = m_message->m_data[kdest];
+			}
+		}
+
+		for (kdest=0; kdest<Kdest.m_K; kdest++)
+		{
+			m_message->m_data[kdest] -= vMin;
+		}
+	}
+	else
+	{
+		assert(0);
+	}
+
+	return vMin;
+}
+
+inline void TypeGeneral::Edge::AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir)
+{
+	assert(ksource>=0 && ksource<Ksource.m_K);
+
+	int k;
+
+	if (m_type == POTTS)
+	{
+		for (k=0; k<ksource; k++)
+		{
+			dest->m_data[k] += ((EdgePotts*)this)->m_lambdaPotts;
+		}
+		for (k++; k<Kdest.m_K; k++)
+		{
+			dest->m_data[k] += ((EdgePotts*)this)->m_lambdaPotts;
+		}
+	}
+	else if (m_type == GENERAL)
+	{
+		REAL* data = ((EdgeGeneral*)this)->m_data;
+
+		if (dir == ((EdgeGeneral*)this)->m_dir)
+		{
+			for (k=0; k<Kdest.m_K; k++)
+			{
+				dest->m_data[k] += data[ksource + k*Ksource.m_K];
+			}
+		}
+		else
+		{
+			for (k=0; k<Kdest.m_K; k++)
+			{
+				dest->m_data[k] += data[k + ksource*Kdest.m_K];
+			}
+		}
+	}
+	else
+	{
+		assert(0);
+	}
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+#endif

libs/Math/TRWS/typePotts.h

@@ -0,0 +1,443 @@
+/******************************************************************
+typePotts.h
+
+Energy function with Potts interactions:
+   E(x)   =   \sum_i D_i(x_i)   +   \sum_ij V_ij(x_i,x_j)
+   where x_i \in {0, 1, ..., K-1},
+   V_ij(ki, kj) = 0 if ki==kj, and lambda_ij otherwise.
+   lambda_ij must be non-negative.
+
+Example usage:
+
+Minimize function E(x,y) = Dx(x) + Dy(y) + lambda*[x != y] where 
+  x,y \in {0,1,2},
+  Dx(0) = 0, Dx(1) = 1, Dx(2) = 2,
+  Dy(0) = 3, Dy(1) = 4, Dy(2) = 5,
+  lambda = 6,
+  [.] is 1 if it's argument is true, and 0 otherwise.
+
+
+
+#include <stdio.h>
+#include "MRFEnergy.h"
+
+void testPotts()
+{
+	MRFEnergy<TypePotts>* mrf;
+	MRFEnergy<TypePotts>::NodeId* nodes;
+	MRFEnergy<TypePotts>::Options options;
+	TypePotts::REAL energy, lowerBound;
+
+	const int nodeNum = 2; // number of nodes
+	const int K = 3; // number of labels
+	TypePotts::REAL D[K];
+	int x, y;
+
+	mrf = new MRFEnergy<TypePotts>(TypePotts::GlobalSize(K));
+	nodes = new MRFEnergy<TypePotts>::NodeId[nodeNum];
+
+	// construct energy
+	D[0] = 0; D[1] = 1; D[2] = 2;
+	nodes[0] = mrf->AddNode(TypePotts::LocalSize(), TypePotts::NodeData(D));
+	D[0] = 3; D[1] = 4; D[2] = 5;
+	nodes[1] = mrf->AddNode(TypePotts::LocalSize(), TypePotts::NodeData(D));
+	mrf->AddEdge(nodes[0], nodes[1], TypePotts::EdgeData(6));
+
+	// 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\n", x, y);
+
+	// done
+	delete nodes;
+	delete mrf;
+}
+
+*******************************************************************/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#ifndef __TYPEPOTTS_H__
+#define __TYPEPOTTS_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypePotts
+{
+private:
+	struct Vector; // node parameters and messages
+	struct Edge; // stores edge information and either forward or backward message
+
+public:
+	// types declarations
+	typedef int Label;
+	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 K);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		int		m_K; // number of labels
+	};
+
+	struct LocalSize // number of labels is stored at MRFEnergy::m_Kglobal
+	{
+	};
+
+	struct NodeData
+	{
+		NodeData(const REAL* data); // data = pointer to array of size MRFEnergy::m_Kglobal
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		const REAL*	m_data;
+	};
+
+	struct EdgeData
+	{
+		EdgeData(REAL lambdaPotts);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_lambdaPotts;
+	};
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	////////////////////////// Visible only to MRFEnergy /////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+
+private:
+friend class MRFEnergy<TypePotts>;
+
+	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:
+		// edge information
+		REAL		m_lambdaPotts;
+
+		// message
+		Vector		m_message;
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+inline TypePotts::GlobalSize::GlobalSize(int K)
+{
+	m_K = K;
+}
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypePotts::NodeData::NodeData(const REAL* data)
+{
+	m_data = data;
+}
+
+inline TypePotts::EdgeData::EdgeData(REAL lambdaPotts)
+{
+	m_lambdaPotts = lambdaPotts;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypePotts::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	if (Kglobal.m_K < 1)
+	{
+		return -1;
+	}
+	return Kglobal.m_K*sizeof(REAL);
+}
+inline void TypePotts::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	memcpy(m_data, data.m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypePotts::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 TypePotts::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	memset(m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypePotts::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	memcpy(m_data, V->m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypePotts::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 TypePotts::REAL TypePotts::Vector::GetValue(GlobalSize Kglobal, LocalSize K, Label k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline TypePotts::REAL TypePotts::Vector::ComputeMin(GlobalSize Kglobal, LocalSize K, Label& kMin)
+{
+	REAL vMin = m_data[0];
+	kMin = 0;
+	for (int k=1; k<Kglobal.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+			kMin = k;
+		}
+	}
+
+	return vMin;
+}
+
+inline TypePotts::REAL TypePotts::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 TypePotts::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return Kglobal.m_K;
+}
+
+inline TypePotts::REAL TypePotts::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline void TypePotts::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 TypePotts::Edge::GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data)
+{
+	if (data.m_lambdaPotts < 0)
+	{
+		return -1;
+	}
+	return sizeof(Edge) - sizeof(Vector) + Kglobal.m_K*sizeof(REAL);
+}
+
+inline int TypePotts::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	return 0;
+}
+
+inline void TypePotts::Edge::Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj)
+{
+	m_lambdaPotts = data.m_lambdaPotts;
+	memset(m_message.m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline TypePotts::Vector* TypePotts::Edge::GetMessagePtr()
+{
+	return &m_message;
+}
+
+inline void TypePotts::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+}
+
+inline TypePotts::REAL TypePotts::Edge::UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* buf)
+{
+	int k;
+	REAL vMin;
+
+	m_message.m_data[0] = gamma*source->m_data[0] - m_message.m_data[0];
+	vMin = m_message.m_data[0];
+
+	for (k=1; k<Kglobal.m_K; k++)
+	{
+		m_message.m_data[k] = gamma*source->m_data[k] - m_message.m_data[k];
+		if (vMin > m_message.m_data[k])
+		{
+			vMin = m_message.m_data[k];
+		}
+	}
+
+	for (k=0; k<Kglobal.m_K; k++)
+	{
+		m_message.m_data[k] -= vMin;
+		if (m_message.m_data[k] > m_lambdaPotts)
+		{
+			m_message.m_data[k] = m_lambdaPotts;
+		}
+	}
+
+	return vMin;
+}
+
+inline void TypePotts::Edge::AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir)
+{
+	assert(ksource>=0 && ksource<Kglobal.m_K);
+
+	int k;
+
+	for (k=0; k<ksource; k++)
+	{
+		dest->m_data[k] += m_lambdaPotts;
+	}
+	for (k++; k<Kglobal.m_K; k++)
+	{
+		dest->m_data[k] += m_lambdaPotts;
+	}
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+#endif

libs/Math/TRWS/typeTruncatedLinear.h

@@ -0,0 +1,463 @@
+/******************************************************************
+typeTruncatedLinear.h
+
+Energy function with truncated linear interactions:
+   E(x)   =   \sum_i D_i(x_i)   +   \sum_ij V_ij(x_i,x_j)
+   where x_i \in {0, 1, ..., K-1},
+   V_ij(ki, kj) = min { alpha_ij*|ki-kj|, lambda_ij }.
+   alpha_ij and lambda_ij must be non-negative.
+
+Example usage:
+
+Minimize function E(x,y) = Dx(x) + Dy(y) + min { alpha*|x - y| , lambda } where 
+  x,y \in {0,1,2},
+  Dx(0) = 0, Dx(1) = 1, Dx(2) = 2,
+  Dy(0) = 3, Dy(1) = 4, Dy(2) = 5,
+  alpha = 6,
+  lambda = 7
+
+
+
+
+#include <stdio.h>
+#include "MRFEnergy.h"
+
+void testTruncatedLinear()
+{
+	MRFEnergy<TypeTruncatedLinear>* mrf;
+	MRFEnergy<TypeTruncatedLinear>::NodeId* nodes;
+	MRFEnergy<TypeTruncatedLinear>::Options options;
+	TypeTruncatedLinear::REAL energy, lowerBound;
+
+	const int nodeNum = 2; // number of nodes
+	const int K = 3; // number of labels
+	TypeTruncatedLinear::REAL D[K];
+	int x, y;
+
+	mrf = new MRFEnergy<TypeTruncatedLinear>(TypeTruncatedLinear::GlobalSize(K));
+	nodes = new MRFEnergy<TypeTruncatedLinear>::NodeId[nodeNum];
+
+	// construct energy
+	D[0] = 0; D[1] = 1; D[2] = 2;
+	nodes[0] = mrf->AddNode(TypeTruncatedLinear::LocalSize(), TypeTruncatedLinear::NodeData(D));
+	D[0] = 3; D[1] = 4; D[2] = 5;
+	nodes[1] = mrf->AddNode(TypeTruncatedLinear::LocalSize(), TypeTruncatedLinear::NodeData(D));
+	mrf->AddEdge(nodes[0], nodes[1], TypeTruncatedLinear::EdgeData(6, 7));
+
+	// 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\n", x, y);
+
+	// done
+	delete nodes;
+	delete mrf;
+}
+
+*******************************************************************/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#ifndef __TYPETRUNCATEDLINEAR_H__
+#define __TYPETRUNCATEDLINEAR_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypeTruncatedLinear
+{
+private:
+	struct Vector; // node parameters and messages
+	struct Edge; // stores edge information and either forward or backward message
+
+public:
+	// types declarations
+	typedef int Label;
+	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 K);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		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 alpha, REAL lambda);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_alpha;
+		REAL		m_lambda;
+	};
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	////////////////////////// Visible only to MRFEnergy /////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+
+private:
+friend class MRFEnergy<TypeTruncatedLinear>;
+
+	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:
+		// edge information
+		REAL		m_alpha;
+		REAL		m_lambda;
+
+		// message
+		Vector		m_message;
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+inline TypeTruncatedLinear::GlobalSize::GlobalSize(int K)
+{
+	m_K = K;
+}
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypeTruncatedLinear::NodeData::NodeData(REAL* data)
+{
+	m_data = data;
+}
+
+inline TypeTruncatedLinear::EdgeData::EdgeData(REAL alpha, REAL lambda)
+{
+	m_alpha = alpha;
+	m_lambda = lambda;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypeTruncatedLinear::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	if (Kglobal.m_K < 1)
+	{
+		return -1;
+	}
+	return Kglobal.m_K*sizeof(REAL);
+}
+inline void TypeTruncatedLinear::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	memcpy(m_data, data.m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedLinear::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 TypeTruncatedLinear::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	memset(m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedLinear::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	memcpy(m_data, V->m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedLinear::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 TypeTruncatedLinear::REAL TypeTruncatedLinear::Vector::GetValue(GlobalSize Kglobal, LocalSize K, Label k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline TypeTruncatedLinear::REAL TypeTruncatedLinear::Vector::ComputeMin(GlobalSize Kglobal, LocalSize K, Label& kMin)
+{
+	REAL vMin = m_data[0];
+	kMin = 0;
+	for (int k=1; k<Kglobal.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+			kMin = k;
+		}
+	}
+
+	return vMin;
+}
+
+inline TypeTruncatedLinear::REAL TypeTruncatedLinear::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 TypeTruncatedLinear::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return Kglobal.m_K;
+}
+
+inline TypeTruncatedLinear::REAL TypeTruncatedLinear::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline void TypeTruncatedLinear::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 TypeTruncatedLinear::Edge::GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data)
+{
+	if (data.m_alpha < 0 || data.m_lambda < 0)
+	{
+		return -1;
+	}
+	return sizeof(Edge) - sizeof(Vector) + Kglobal.m_K*sizeof(REAL);
+}
+
+inline int TypeTruncatedLinear::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	return 0;
+}
+
+inline void TypeTruncatedLinear::Edge::Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj)
+{
+	m_alpha = data.m_alpha;
+	m_lambda = data.m_lambda;
+	memset(m_message.m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline TypeTruncatedLinear::Vector* TypeTruncatedLinear::Edge::GetMessagePtr()
+{
+	return &m_message;
+}
+
+inline void TypeTruncatedLinear::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+}
+
+inline TypeTruncatedLinear::REAL TypeTruncatedLinear::Edge::UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* buf)
+{
+	int k;
+	REAL vMin;
+
+	m_message.m_data[0] = gamma*source->m_data[0] - m_message.m_data[0];
+	vMin = m_message.m_data[0];
+
+	for (k=1; k<Kglobal.m_K; k++)
+	{
+		m_message.m_data[k] = gamma*source->m_data[k] - m_message.m_data[k];
+		if (m_message.m_data[k] > m_message.m_data[k-1] + m_alpha)
+		{
+			m_message.m_data[k] = m_message.m_data[k-1] + m_alpha;
+		}
+		else if (vMin > m_message.m_data[k])
+		{
+			vMin = m_message.m_data[k];
+		}
+	}
+
+	k--;
+	m_message.m_data[k] -= vMin;
+	if (m_message.m_data[k] > m_lambda)
+	{
+		m_message.m_data[k] = m_lambda;
+	}
+
+	for (k--; k>=0; k--)
+	{
+		m_message.m_data[k] -= vMin;
+		if (m_message.m_data[k] > m_message.m_data[k+1] + m_alpha)
+		{
+			m_message.m_data[k] = m_message.m_data[k+1] + m_alpha;
+		}
+		if (m_message.m_data[k] > m_lambda)
+		{
+			m_message.m_data[k] = m_lambda;
+		}
+	}
+
+	return vMin;
+}
+
+inline void TypeTruncatedLinear::Edge::AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir)
+{
+	assert(ksource>=0 && ksource<Kglobal.m_K);
+
+	int k;
+
+	for (k=0; k<ksource; k++)
+	{
+		dest->m_data[k] += (ksource-k)*m_alpha < m_lambda ? (ksource-k)*m_alpha : m_lambda;
+	}
+	for (k++; k<Kglobal.m_K; k++)
+	{
+		dest->m_data[k] += m_alpha*(k-ksource) < m_lambda ? m_alpha*(k-ksource) : m_lambda;
+	}
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+#endif

libs/Math/TRWS/typeTruncatedLinear2D.h

@@ -0,0 +1,491 @@
+/******************************************************************
+typeTruncatedLinear2D.h
+
+Energy function with 2-dimensional truncated linear 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)| + alpha2_ij*|ki(2)-kj(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)| + alpha2*|x(2) - y(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 testTruncatedLinear2D()
+{
+	MRFEnergy<TypeTruncatedLinear2D>* mrf;
+	MRFEnergy<TypeTruncatedLinear2D>::NodeId* nodes;
+	MRFEnergy<TypeTruncatedLinear2D>::Options options;
+	TypeTruncatedLinear2D::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
+	TypeTruncatedLinear2D::REAL D[K];
+	TypeTruncatedLinear2D::Label x, y;
+
+	mrf = new MRFEnergy<TypeTruncatedLinear2D>(TypeTruncatedLinear2D::GlobalSize(KX, KY));
+	nodes = new MRFEnergy<TypeTruncatedLinear2D>::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(TypeTruncatedLinear2D::LocalSize(), TypeTruncatedLinear2D::NodeData(D));
+	D[0] = 0; D[1] = 0; D[2] = 0;
+	D[3] = 0; D[4] = 0; D[5] = 0;
+	nodes[1] = mrf->AddNode(TypeTruncatedLinear2D::LocalSize(), TypeTruncatedLinear2D::NodeData(D));
+	mrf->AddEdge(nodes[0], nodes[1], TypeTruncatedLinear2D::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 __TYPETRUNCATEDLINEAR2D_H__
+#define __TYPETRUNCATEDLINEAR2D_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypeTruncatedLinear2D
+{
+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<TypeTruncatedLinear2D>;
+
+	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:
+		// edge information
+		REAL		m_alphaX, m_alphaY;
+		REAL		m_lambda;
+
+		// message
+		Vector		m_message;
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+inline TypeTruncatedLinear2D::GlobalSize::GlobalSize(int KX, int KY)
+{
+	m_KX = KX;
+	m_KY = KY;
+	m_K = KX*KY;
+}
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypeTruncatedLinear2D::NodeData::NodeData(REAL* data)
+{
+	m_data = data;
+}
+
+inline TypeTruncatedLinear2D::EdgeData::EdgeData(REAL alphaX, REAL alphaY, REAL lambda)
+{
+	m_alphaX = alphaX;
+	m_alphaY = alphaY;
+	m_lambda = lambda;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypeTruncatedLinear2D::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	if (Kglobal.m_KX < 1 || Kglobal.m_KY < 2)
+	{
+		return -1;
+	}
+	return Kglobal.m_K*sizeof(REAL);
+}
+inline void TypeTruncatedLinear2D::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	memcpy(m_data, data.m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedLinear2D::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 TypeTruncatedLinear2D::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	memset(m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedLinear2D::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	memcpy(m_data, V->m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedLinear2D::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 TypeTruncatedLinear2D::REAL TypeTruncatedLinear2D::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 TypeTruncatedLinear2D::REAL TypeTruncatedLinear2D::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 TypeTruncatedLinear2D::REAL TypeTruncatedLinear2D::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 TypeTruncatedLinear2D::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return Kglobal.m_K;
+}
+
+inline TypeTruncatedLinear2D::REAL TypeTruncatedLinear2D::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline void TypeTruncatedLinear2D::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 TypeTruncatedLinear2D::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 TypeTruncatedLinear2D::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	return 0;
+}
+
+inline void TypeTruncatedLinear2D::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 TypeTruncatedLinear2D::Vector* TypeTruncatedLinear2D::Edge::GetMessagePtr()
+{
+	return &m_message;
+}
+
+inline void TypeTruncatedLinear2D::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+}
+
+inline TypeTruncatedLinear2D::REAL TypeTruncatedLinear2D::Edge::UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* buf)
+{
+	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 (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++, sourcePtr++, destPtr++)
+	{
+		destPtr[0] = gamma*sourcePtr[0] - destPtr[0];
+		if (k.m_kx > 0 && destPtr[0] > destPtr[-1] + m_alphaX)
+		{
+			destPtr[0] = destPtr[-1] + m_alphaX;
+		}
+		if (k.m_ky > 0 && destPtr[0] > destPtr[-Kglobal.m_KX] + m_alphaY)
+		{
+			destPtr[0] = destPtr[-Kglobal.m_KX] + m_alphaY;
+		}
+		else if (vMin > destPtr[0])
+		{
+			vMin = destPtr[0];
+		}
+	}
+
+	destPtr--;
+	for (k.m_ky=Kglobal.m_KY-1; k.m_ky>=0; k.m_ky--)
+	for (k.m_kx=Kglobal.m_KX-1; k.m_kx>=0; k.m_kx--, destPtr--)
+	{
+		destPtr[0] -= vMin;
+		if (k.m_kx < Kglobal.m_KX-1 && destPtr[0] > destPtr[1] + m_alphaX)
+		{
+			destPtr[0] = destPtr[1] + m_alphaX;
+		}
+		if (k.m_ky < Kglobal.m_KY-1 && destPtr[0] > destPtr[Kglobal.m_KX] + m_alphaY)
+		{
+			destPtr[0] = destPtr[Kglobal.m_KX] + m_alphaY;
+		}
+		if (destPtr[0] > m_lambda)
+		{
+			destPtr[0] = m_lambda;
+		}
+	}
+
+	return vMin;
+}
+
+inline void TypeTruncatedLinear2D::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 : ksource.m_kx - k.m_kx)
+		          + m_alphaY*((k.m_ky > ksource.m_ky) ? k.m_ky - ksource.m_ky : ksource.m_ky - k.m_ky);
+		destPtr[0] += (cost < m_lambda) ? cost : m_lambda;
+	}
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+#endif

libs/Math/TRWS/typeTruncatedQuadratic.h

@@ -0,0 +1,533 @@
+/******************************************************************
+typeTruncatedQuadratic.h
+
+Energy function with 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, ..., K-1},
+   V_ij(ki, kj) = min { alpha_ij*(ki-kj)^2, lambda_ij }.
+   alpha_ij and lambda_ij must be non-negative.
+
+Example usage:
+
+Minimize function E(x,y) = Dx(x) + Dy(y) + min { alpha*(x - y)^2 , lambda } where 
+  x,y \in {0,1,2},
+  Dx(0) = 0, Dx(1) = 1, Dx(2) = 2,
+  Dy(0) = 3, Dy(1) = 4, Dy(2) = 5,
+  alpha = 6,
+  lambda = 7
+
+
+
+
+#include <stdio.h>
+#include "MRFEnergy.h"
+
+void testTruncatedQuadratic()
+{
+	MRFEnergy<TypeTruncatedQuadratic>* mrf;
+	MRFEnergy<TypeTruncatedQuadratic>::NodeId* nodes;
+	MRFEnergy<TypeTruncatedQuadratic>::Options options;
+	TypeTruncatedQuadratic::REAL energy, lowerBound;
+
+	const int nodeNum = 2; // number of nodes
+	const int K = 3; // number of labels
+	TypeTruncatedQuadratic::REAL D[K];
+	int x, y;
+
+	mrf = new MRFEnergy<TypeTruncatedQuadratic>(TypeTruncatedQuadratic::GlobalSize(K));
+	nodes = new MRFEnergy<TypeTruncatedQuadratic>::NodeId[nodeNum];
+
+	// construct energy
+	D[0] = 0; D[1] = 1; D[2] = 2;
+	nodes[0] = mrf->AddNode(TypeTruncatedQuadratic::LocalSize(), TypeTruncatedQuadratic::NodeData(D));
+	D[0] = 3; D[1] = 4; D[2] = 5;
+	nodes[1] = mrf->AddNode(TypeTruncatedQuadratic::LocalSize(), TypeTruncatedQuadratic::NodeData(D));
+	mrf->AddEdge(nodes[0], nodes[1], TypeTruncatedQuadratic::EdgeData(6, 7));
+
+	// 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\n", x, y);
+
+	// done
+	delete nodes;
+	delete mrf;
+}
+
+*******************************************************************/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#ifndef __TYPETRUNCATEDQUADRATIC_H__
+#define __TYPETRUNCATEDQUADRATIC_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypeTruncatedQuadratic
+{
+private:
+	struct Vector; // node parameters and messages
+	struct Edge; // stores edge information and either forward or backward message
+
+public:
+	// types declarations
+	typedef int Label;
+	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 K);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		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 alpha, REAL lambda);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_alpha;
+		REAL		m_lambda;
+	};
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	////////////////////////// Visible only to MRFEnergy /////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+
+private:
+friend class MRFEnergy<TypeTruncatedQuadratic>;
+
+	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, REAL* source, REAL* dest, int* parabolas, int* intersections);
+
+		// edge information
+		REAL		m_alpha;
+		REAL		m_lambda;
+
+		// message
+		Vector		m_message;
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+inline TypeTruncatedQuadratic::GlobalSize::GlobalSize(int K)
+{
+	m_K = K;
+}
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypeTruncatedQuadratic::NodeData::NodeData(REAL* data)
+{
+	m_data = data;
+}
+
+inline TypeTruncatedQuadratic::EdgeData::EdgeData(REAL alpha, REAL lambda)
+{
+	m_alpha = alpha;
+	m_lambda = lambda;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypeTruncatedQuadratic::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	if (Kglobal.m_K < 1)
+	{
+		return -1;
+	}
+	return Kglobal.m_K*sizeof(REAL);
+}
+inline void TypeTruncatedQuadratic::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	memcpy(m_data, data.m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic::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 TypeTruncatedQuadratic::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	memset(m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	memcpy(m_data, V->m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic::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 TypeTruncatedQuadratic::REAL TypeTruncatedQuadratic::Vector::GetValue(GlobalSize Kglobal, LocalSize K, Label k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline TypeTruncatedQuadratic::REAL TypeTruncatedQuadratic::Vector::ComputeMin(GlobalSize Kglobal, LocalSize K, Label& kMin)
+{
+	REAL vMin = m_data[0];
+	kMin = 0;
+	for (int k=1; k<Kglobal.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+			kMin = k;
+		}
+	}
+
+	return vMin;
+}
+
+inline TypeTruncatedQuadratic::REAL TypeTruncatedQuadratic::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 TypeTruncatedQuadratic::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return Kglobal.m_K;
+}
+
+inline TypeTruncatedQuadratic::REAL TypeTruncatedQuadratic::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline void TypeTruncatedQuadratic::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 TypeTruncatedQuadratic::Edge::GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data)
+{
+	if (data.m_alpha < 0 || data.m_lambda < 0)
+	{
+		return -1;
+	}
+	return sizeof(Edge) - sizeof(Vector) + Kglobal.m_K*sizeof(REAL);
+}
+
+inline int TypeTruncatedQuadratic::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	int K = vectorMaxSizeInBytes / sizeof(REAL);
+	return K*sizeof(REAL) + (2*K+1)*sizeof(int);
+}
+
+inline void TypeTruncatedQuadratic::Edge::Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj)
+{
+	m_alpha = data.m_alpha;
+	m_lambda = data.m_lambda;
+	memset(m_message.m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline TypeTruncatedQuadratic::Vector* TypeTruncatedQuadratic::Edge::GetMessagePtr()
+{
+	return &m_message;
+}
+
+inline void TypeTruncatedQuadratic::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+}
+
+inline TypeTruncatedQuadratic::REAL TypeTruncatedQuadratic::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)));
+	int k;
+	REAL vMin;
+
+	buf[0] = gamma*source->m_data[0] - m_message.m_data[0];
+	vMin = buf[0];
+
+	for (k=1; k<Kglobal.m_K; k++)
+	{
+		buf[k] = gamma*source->m_data[k] - m_message.m_data[k];
+		if (vMin > buf[k])
+		{
+			vMin = buf[k];
+		}
+	}
+
+	if (m_alpha == 0)
+	{
+		for (k=0; k<Kglobal.m_K; k++)
+		{
+			m_message.m_data[k] = 0;
+		}
+	}
+	else
+	{
+		DistanceTransformL2(Kglobal.m_K, buf, m_message.m_data, parabolas, intersections);
+
+		for (k=0; k<Kglobal.m_K; k++)
+		{
+			m_message.m_data[k] -= vMin;
+			if (m_message.m_data[k] > m_lambda)
+			{
+				m_message.m_data[k] = m_lambda;
+			}
+		}
+	}
+
+	return vMin;
+}
+
+inline void TypeTruncatedQuadratic::Edge::AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir)
+{
+	assert(ksource>=0 && ksource<Kglobal.m_K);
+
+	int k;
+
+	for (k=0; k<Kglobal.m_K; k++)
+	{
+		dest->m_data[k] += m_alpha*(k-ksource)*(k-ksource) < m_lambda ? m_alpha*(k-ksource)*(k-ksource) : m_lambda;
+	}
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+
+
+
+inline void TypeTruncatedQuadratic::Edge::DistanceTransformL2(int K, REAL* source, REAL* dest, int* parabolas, int* intersections)
+{
+	assert(m_alpha > 0);
+
+	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[j] - source[i]) / (m_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[i] = source[j] + m_alpha*(i-j)*(i-j);
+		}
+	}
+}
+
+#endif

libs/Math/TRWS/typeTruncatedQuadratic2D.h

@@ -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