function [belief2d, pred_labels] = bp_infer(features, edgemat, theta, tau, num_bases, bp_pars)

w = theta(1:num_bases);
b = theta(num_bases+1);
v = theta(num_bases+2:end-1);
c = theta(end);

[dim_feature, num_rows, num_cols] = size(features);

node_features = reshape(features, dim_feature, num_rows*num_cols);
node_potentials_p = 1 ./ (1 + exp(-w' * node_features - b));
node_potentials_n = 1 ./ (1 + exp(w' * node_features + b));
node_potentials = [node_potentials_p' node_potentials_n'];

num_edges = size(edgemat,1);
edge_features = zeros(dim_feature, num_edges);
for i = 1 : size(edgemat,1)
    r1 = edgemat(i, 1); c1 = edgemat(i, 2); 
    r2 = edgemat(i, 3); c2 = edgemat(i, 4);
    edge_features(:, i) = features(:, r1, c1) - features(:, r2, c2);
end

edge_potentials_pp = 1 ./ (1 + exp(-v' * (edge_features .* edge_features) - c));
edge_potentials_pn = 1 ./ (1 + exp(v' * (edge_features .* edge_features) + c));
edge_potentials_np = edge_potentials_pn;
edge_potentials_nn = edge_potentials_pp;

edge_potentials = zeros(2, 2, num_edges);
edge_potentials(1, 1, :) = edge_potentials_pp;
edge_potentials(1, 2, :) = edge_potentials_pn;
edge_potentials(2, 1, :) = edge_potentials_np;
edge_potentials(2, 2, :) = edge_potentials_nn;

node_potentials = double(node_potentials);
edge_potentials = double(edge_potentials.^tau);

[V, E] = StarEdge_MakeEdgeNums_Lattice2(num_rows, num_cols);

[node_belief, num_iter, msgs, edge_belief] = BP_General_C(...
        edge_potentials, node_potentials, bp_pars.maximize, bp_pars.max_iters,...
        bp_pars.tol, int32(V), int32(E));
    
belief2d = reshape(node_belief, [num_rows num_cols 2]);

[junk, pred_labels] = max(belief2d, [], 3);
pred_labels = -pred_labels*2 + 3;
    
    
    
    
    
    
    
    
    
    