	clear; 
	clc; 
	close all;

	% Pairwise potential parameters
	%
	lambda = 20;
	disc_K = 2;

	% FastPD params
	%
	max_iters = 30;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	% Compute unary potentials
	%
	im = imread('Tsukuba_left.ppm');
	sizex = size(im,2);
	sizey = size(im,1);
	numlabels = 16;

	fp = fopen( 'tsukuba_unary.bin', 'rb' );
	unary = fread( fp, inf, 'int32' );
	unary = reshape( unary, sizex*sizey, numlabels );
	fclose(fp);
	
	% Generate the horizontal edges of the MRF graph (in this case MRF-graph = 4-connected image grid)
	%
	[x0,y0] = meshgrid(1:sizex-1,1:sizey);
	hx = x0+1; hy = y0;
	hor_inds0 = sub2ind( [sizey sizex], y0(:), x0(:) )' - 1; % I assume C-indexing 
	hor_inds1 = sub2ind( [sizey sizex], hy(:), hx(:) )' - 1;

	% Generate the vertical edges of the MRF graph (in this case MRF-graph = 4-connected image grid) 
	%
	[x1,y1] = meshgrid(1:sizex,1:sizey-1);
	vx = x1; vy = y1+1;
	ver_inds0 = sub2ind( [sizey sizex], y1(:), x1(:) )' - 1;
	ver_inds1 = sub2ind( [sizey sizex], vy(:), vx(:) )' - 1;

	% Concatenate horizontal and vertical edges
	%
	edges = [ hor_inds0 ver_inds0 ;
	          hor_inds1 ver_inds1 ];

	% Set the edge weights
	%
	num_edges = size(edges,2);
	edge_weights = lambda*ones(num_edges,1);

	% Set the pairwise distance (pairwise_potential = edge_weight*pairwise_distance)
	%
	[l0,l1] = meshgrid( 0:numlabels-1, 0:numlabels-1 );
	dist = min( abs(l0-l1), disc_K*ones(numlabels,numlabels) );

	% Run the algorithm
	%
	unary = round(unary);
	edge_weights = round(edge_weights);
	dist = round(dist);

	labeling = FastPD_wrapper( unary, edges, edge_weights, dist, max_iters );
	labeling = reshape(labeling,sizey,sizex);
	figure; imshow(labeling,[]);

