function [ef,s,tmbs] = energy_full(output_image,input_image,lambda)

% Input and output images must be same size

   [oSy,oSx] = size(output_image);
   [iSy,iSx] = size(input_image);

   %if (oSy ~= iSx) || (oSx ~= iSx)
   if (oSy ~= iSy) || (oSx ~= iSx)
    
       error('sizes of output and input images do not match')

   end

   N = oSy;
   M = oSx;

    output_image=double(output_image);
    input_image=double(input_image);

%%%%% set up for computing the TV term %%%%

  a = output_image;

 % pad the image

   b(3:N+2,3:M+2) = a;

   b(3:N+2,2) = a(:,1);
   b(3:N+2,1) = a(:,2);

   b(3:N+2,M+3) = a(:,M);
   b(3:N+2,M+4) = a(:,M-1);

   b(2,:) = b(3,:);
   b(1,:) = b(4,:);

   b(N+3,:) = b(N+2,:);
   b(N+4,:) = b(N+1,:);

 % initialize the neighborhood sums

   n4 = 0;
   n8 = 0;
   n16 = 0;

 % compute the neighborhood sums

   n4 = n4 + sum(sum(max(a - b(3-1:N+2-1,  3:M+2  ),0)));
   n4 = n4 + sum(sum(max(a - b(3+1:N+2+1,  3:M+2  ),0)));
   n4 = n4 + sum(sum(max(a - b(  3:N+2  ,3-1:M+2-1),0)));
   n4 = n4 + sum(sum(max(a - b(  3:N+2  ,3+1:M+2+1),0)));
 
   n8 = n8 + sum(sum(max(a - b(3-1:N+2-1,3-1:M+2-1),0)));
   n8 = n8 + sum(sum(max(a - b(3+1:N+2+1,3-1:M+2-1),0)));
   n8 = n8 + sum(sum(max(a - b(3-1:N+2-1,3+1:M+2+1),0)));
   n8 = n8 + sum(sum(max(a - b(3+1:N+2+1,3+1:M+2+1),0))); 

   n16 = n16 + sum(sum(max(a - b(3-1:N+2-1,3-2:M+2-2),0)));
   n16 = n16 + sum(sum(max(a - b(3+1:N+2+1,3-2:M+2-2),0)));
   n16 = n16 + sum(sum(max(a - b(3-1:N+2-1,3+2:M+2+2),0)));
   n16 = n16 + sum(sum(max(a - b(3+1:N+2+1,3+2:M+2+2),0)));
   n16 = n16 + sum(sum(max(a - b(3-2:N+2-2,3-1:M+2-1),0)));
   n16 = n16 + sum(sum(max(a - b(3+2:N+2+2,3-1:M+2-1),0)));
   n16 = n16 + sum(sum(max(a - b(3-2:N+2-2,3+1:M+2+1),0)));
   n16 = n16 + sum(sum(max(a - b(3+2:N+2+2,3+1:M+2+1),0)));

% sum up the weighted factors to get the TV (perimeter) term

  w1 = 0.2442;
  w2 = 0.0952;
  w3 = 0.0908;

  TV = w1*n4 + w2*n8 + w3*n16;



% Compute the data fidelity (symmetric difference) term

  sym_dif = max(output_image - input_image,0)+ max(input_image - output_image,0);

  data_fidelity = sum(sum(sym_dif));


% Compute energy

  
     ef = TV + lambda*data_fidelity;
   tmbs = TV;
      s = data_fidelity;