%close all;
%clear all;

f = imread('cameraman.tif');
[m,n]=size(f);
dt = 0.1; %timestepping

fn = imnoise(f,'gaussian',0,0.004);
fn = double(fn); f = double(f);

figure(1)
subplot(1,3,1); image(fn); colormap(gray(256));
title('original image')
subplot(1,3,2); image(fn); colormap(gray(256));
title('total variation')
subplot(1,3,3); image(fn); colormap(gray(256));
title('perona malik')

pause
gamma = 1;
deltar=0.01/3;
fpd = explicitheat_nb(5,1e-6,fn); %smooth noisy image
fxf = forwardx_nb(fpd);
fyf = forwardy_nb(fpd);
absgradf = sqrt(fxf.^2 + fyf.^2)+deltar;
gf = 1./(1+ gamma*absgradf.^2);

deltar = 0.01/3; %regularization for singular denom's
lambda = 1/12.5; alpha = 1e0; u = fn; u2 = u;


for i=1:50000
   uold=u;
   fxu = forwardx_nb(u);
   fyu = forwardy_nb(u);
   absgradu = sqrt(fxu.^2 + alpha*(fyu.^2))+deltar;
    
   gamma = 5;
   deltar=0.01/3;
   upd = explicitheat_nb(5,1e-6,u2); %smooth noisy image
   fxu2s = forwardx_nb(upd);
   fyu2s = forwardy_nb(upd);
   absgradu2s = sqrt(fxu2s.^2 + fyu2s.^2)+deltar;
   gu2s = 1./(1+ gamma*absgradu2s.^2);
    
   fxu2 = forwardx_nb(u2);
   fyu2 = forwardy_nb(u2);
   absgradu2 = sqrt(fxu2.^2 + alpha*(fyu2.^2))+deltar;
   
   %tv
   u = u + dt*( lambda*(fn-u) + ( backwardx_nb(fxu./sqrt(fxu.^2+deltar)) + backwardy_nb(alpha*fyu./sqrt(fyu.^2+deltar))   ) );
   %pm
   
   u2 = u2 + dt*( 20*( backwardx_nb(gu2s.*fxu2./sqrt(fxu2.^2+deltar)) + backwardy_nb(gu2s.*fyu2./sqrt(fyu2.^2+deltar))   ) );
    
   norm(u-uold,'fro');

   i
   if rem(i,25)==0
   subplot(1,3,1); image(fn); colormap(gray(256));
   title('original image')
   subplot(1,3,2); image(u); colormap(gray(256));
   title('total variation')
   subplot(1,3,3); image(u2); colormap(gray(256));
   title('perona malik')
   pause
   end
   
end