function E = energy(k,s,y,mu,sigma,x,beta)


[row_max,col_max] = size(x);

row = s(1);
col = s(2);

E = (1/(2*sigma(k)^2))*(y(row,col)-mu(k))^2;


rowp = row+1;
if rowp>row_max
 rowp=1;
end

rowm = row-1;
if rowm < 1
 rowm=row_max;
end

colp = col+1;
if colp>col_max
 colp=1;
end

colm = col-1;
if colm < 1
 colm=col_max;
end


E = E + beta*( (k~=x(rowp,col)) + (k~=x(rowm,col)) + (k~=x(row,colp)) + (k~=x(row,colm)) );


