function U = potentialEnergy(I)

% Subfunction of the Simulated Annealing program.
% Computes the potential energy of an image using a log (laplacian of
% gaussian) filter with some enhancement.


% Initialize auxiliary variables
[nx, ny] = size(I);

% Extend I: Add symmetric borders in order to avoid effects of the borders
% when applying the log filter
border = 10;
extI = zeros(nx+2*border,ny+2*border);
extI((border+1):(nx+border),(border+1):(ny+border)) = I;
for i = 1:border,
    extI(i,:) = [ zeros(1,border) I(border + 2 - i,:) zeros(1,border) ];
    extI(end + 1 - i,:) = [ zeros(1,border) I(end - border - 1 + i,:) zeros(1,border) ];
end
for i = 1:border,
    extI(:,i) = extI(:, 2 * ( border + 1 ) - i );
    extI(:,end + 1 - i) = extI(:,end - 2 * border - 1 + i );
end

% Design and apply the log (laplacian of gaussian) filter
h = fspecial('log', 25, 3.0);
U = filter2(h,extI);

% Cut out the original part of the filtered image
U = U((border + 1):end - border, (border + 1):end - border);

% Cut off the upper 20% (saturate the image) 
N = nx * ny;
u = reshape(U,N,1);
v = sort(u);
cutoff = v(round(.8*N));

for i = 1:N,
    if ( u(i) > cutoff )
        u(i) = cutoff;
    end
end

% Rescale U so that it takes values between 0 and 1
U = reshape(u,nx,ny);
U = U - min(min(U));
U = U ./ max(max(U));

% Apply a non-linear scaling
u = reshape(U,N,1);
for i = 1:N,
    u(i) = u(i)^6;
end

% Return the filtered image
U = reshape(u,nx,ny);