%%% Shock filter brain image example
%%% By Guy Gilboa
%%% 3 shock filters are compared:
%%% Osher-Rudin (original)[1], Alvarez-Mazorra [2], Gilboa et al. [3]

%% Degradatin of original image
ro=0.5;  % Gaussian blurring kernel std
std_n=10;  % std of white Gaussian noise

[I,map]=imread('shock_brain_orig.bmp');
I=ind2gray(I,map)*256;
figure(1); imshow(uint8(I))

%Blur
Ig = gauss(I,11,ro^2);

% Add noise
%std_n=0; % no noise
randn('seed',1259);
noise=randn(size(I))*std_n;
Ign=Ig+noise;
snr=db(Ig,Ign)
figure(2); clf; imshow(uint8(Ign))

%% time and spatial steps
dt=0.1; h=1;

% Original shock filter [1]
iter=10;
J=shock(Ign,iter,dt,h,'org');  
figure(3); clf;imshow(uint8(J))

% Alvarez-Mazorra regularization [2] 
% parameters
iter=20;  % iterations
sig2=4;  % sigma^2 of Gaussian conv. of second derivative 
C=0.2;  % diffusion in level-sets direction

% evolution
J=shock(Ign,iter,dt,h,'alv',[C,sig2]);   
figure(4); clf;imshow(uint8(J))

% Complex shock filter [3]
% parameters
iter=30;      % iterations
lam=0.1; % complex diffusion across the edge 
lam_tld=0.2;      % real diffusion in level-sets direction
a = 2;  % slope of "soft sign"
theta=pi/1000;  % phase angle of complex part (default value)
% evolution
J=shock(Ign,iter,dt,h,'cmp',[lam,lam_tld,a]);   
figure(5); imshow(uint8(real(J)))
figure(6); imshow(uint8(4*imag(J)/theta+128))

% Refs.
% [1] S.J. Osher and L. I. Rudin, "Feature-Oriented Image enhancement using Shock Filters", SIAM J. Numer. Anal. 27, pp. 919-940, 1990.
% [2] L. Alvarez and L. Mazorra. "Signal and image restoration using shock filters and anisotropic diffusion". SIAM J. Numer. Anal., 31(2):590-605, 1994.
% [3] G. Gilboa, N. Sochen, Y.Y. Zeevi, "Regularized shock filters and complex diffusion",  - ECCV-'02, LNCS 2350, pp. 399-313, Springer-Verlag 2002.