function nonlinear(image1);
% performs edge detection using artificial dissipation
% Notice that my output images are not binary as in the article.

%clear

%image=imread('barbara.jpg');  

image2=double(image1);
image =full(image2);

epsilon2=0.5;
epsilon4=1/32;

D2x = [0 0 0;1 -2 1;0 0 0];
D4x = [0 0 0 0 0;0 0 0 0 0;1 -4 6 -4 1;0 0 0 0 0;0 0 0 0 0];
D2y = [0 1 0;0 -2 0;0 1 0];
D4y = [0 0 1 0 0;0 0 -4 0 0;0 0 6 0 0;0 0 -4 0 0;0 0 1 0 0];
D2m = [0 0 1;0 -2 0;1 0 0];
D4m = [0 0 0 0 1;0 0 0 -4 0;0 0 6 0 0;0 -4 0 0 0;1 0 0 0 0];
D2n = [1 0 0;0 -2 0;0 0 1];
D4n = [1 0 0 0 0;0 -4 0 0 0;0 0 6 0 0;0 0 0 -4 0;0 0 0 0 1];

f2x = filter2(D2x,image);
f4x = filter2(D4x,image);
f2y = filter2(D2y,image);
f4y = filter2(D4y,image);
f2m = filter2(D2m,image);
f4m = filter2(D4m,image);
f2n = filter2(D2n,image);
f4n = filter2(D4n,image);

%dx =image(:,[2:end])-image(:,[1:end-1]);
%dy =image([2:end],:)-image([1:end-1],:);
%dx2=dx(:,[2:end])-dx(:,[1:end-1]);
%dy2=dy([2:end],:)-dy([1:end-1],:);
Dpx = [0 0 0;0 -1 1;0 0 0];
Dmx = [0 0 0;-1 1 0;0 0 0];
Dax = [0 0 0;1 2 1;0 0 0];
Dpy = [0 0 0;0 -1 0;0 1 0];
Dmy = [0 -1 0;0 1 0;0 0 0];
Day = [0 1 0;0 2 0;0 1 0];
Dpm = [0 0 0;0 -1 0;1 0 0];
Dmm = [0 0 -1;0 1 0;0 0 0];
Dam = [0 0 1;0 2 0;1 0 0];
Dpn = [0 0 0;0 -1 0;0 0 1];
Dmn = [-1 0 0;0 1 0;0 0 0];
Dan = [1 0 0;0 2 0;0 0 1];
%
fpx= filter2(Dpx,image);
fmx= filter2(Dmx,image);
fax= filter2(Dax,image);
fpy= filter2(Dpy,image);
fmy= filter2(Dmy,image);
fay= filter2(Day,image);
fpm= filter2(Dpm,image);
fmm= filter2(Dmm,image);
fam= filter2(Dam,image);
fpn= filter2(Dpn,image);
fmn= filter2(Dmn,image);
fan= filter2(Dan,image);
%
numx = abs(f2x);
numy = abs(f2y);
numm = abs(f2m);
numn = abs(f2n);
denx = abs(fpx)+ abs(fmx);
deny = abs(fpy)+ abs(fmy);
denm = abs(fpm)+ abs(fmm);
denn = abs(fpn)+ abs(fmn);
enteps =0.5;
epsilon24=2*epsilon2;
epsilons=1.e-5;
denx = enteps*fax + (1-enteps)*denx + epsilons;
deny = enteps*fay + (1-enteps)*deny + epsilons;
denm = enteps*fam + (1-enteps)*denm + epsilons;
denn = enteps*fan + (1-enteps)*denn + epsilons;
dpx = numx./denx;
dpy = numy./deny;
dpm = numm./denm;
dpn = numn./denn;
dppx=dpx(:,[2:end]);
dppy=dpy([2:end],:);
dppm=dpx(:,[2:end]);
dppn=dpy([2:end],:);
[daxx,daxy]=size(dppx);
dppx=[dppx,zeros(daxx,1)];
[dayx,dayy]=size(dppy);
dppy=[zeros(1,dayy); dppy];
[damx,damy]=size(dppm);
dppm=[dppm,zeros(damx,1)];
[danx,dany]=size(dppy);
dppn=[zeros(1,dany); dppn];
epx = max(dpx,dppx);
epy = max(dpy,dppy);
epm = max(dpm,dppm);
epn = max(dpn,dppn);
epm=dpm;
epn=dpn;
e4x = max(epsilon4-epsilon24*epx,0);
e4y = max(epsilon4-epsilon24*epy,0);
e4m = max(epsilon4-epsilon24*epm,0);
e4n = max(epsilon4-epsilon24*epn,0);
%
%dx =image(:,[2:end])-image(:,[1:end-1]);
%dy =image([2:end],:)-image([1:end-1],:);
%dx2=dx(:,[2:end])-dx(:,[1:end-1]);
%dy2=dy([2:end],:)-dy([1:end-1],:);
%
%pr=epx/epsilon2
epxl=1;
epyl=1;
epml=1;
epnl=1;
%pr=emx/epsilon4
e4xl=epsilon4;
e4yl=epsilon4;
e4ml=epsilon4;
e4nl=epsilon4;
%e4x=epsilon4;
%e4y=epsilon4;
%e4m=epsilon4;
%e4n=epsilon4;
%
epx=1;
epy=1;
epm=1;
epn=1;
%
e4x=epsilon4;
e4y=epsilon4;
e4m=epsilon4;
e4n=epsilon4;
%
tmp1 = epsilon2*(epxl*f2x+epyl*f2y+epml*f2m+epnl*f2n) ;
tmp2 = epsilon2*(epx.*f2x+epy.*f2y+epm.*f2m+epn.*f2n) ;
tmp1 = tmp1 - (e4xl*f4x+e4yl*f4y+e4ml*f4m+e4nl*f4n);   
tmp2 = tmp2 - (e4x.*f4x+e4y.*f4y+e4m.*f4m+e4n.*f4n);    

c = canny(image,1);   % see canny.m

new_img  = c + tmp1;
new_img2= c + tmp2;

figure(3) 
imshow(uint8(c))         % the output of canny
title('canny ');
%figure, imshow(uint8(tmp))      % the "correction" 
figure(4)
imshow(uint8(new_img))   % the new image ( canny+correction )
title('correction - linear');
figure(5)
imshow(uint8(new_img2))   % the new image ( canny+correction )
title('correction - nonlinear');