hbtfilt(x)

function hbtfilt(x)
%Reading Image
Im=imread(x);
Im=im2double(Im);
[row1 col1]=size(Im);
%Defining Window
W=2;
Z=1;
ws=2*W+1;

%Extracting the neighborhoods
Im=padarray(Im,[W W],0,'both');
Ng = im2col(Im,[ws ws],'sliding');
Ngc=Ng';
pca=princomp(Ngc);
%Calculating the PCA
e2v=pca(:,2);
e3v=pca(:,3);
e2m=vec2mat(e2v,ws)';
e3m=vec2mat(e3v,ws)';

%Drawing the Mesh
x=-W:1:W;
y=-W:1:W;
figure,mesh(x,y,e2m);
grid off;
figure,mesh(x,y,e3m);
grid off;

%Calculating sigma value
k=1;
for j=0.01:0.01:5
  sig(k)=j;
  h1m=hbtfilter(j,W);
  h2m=rot90(h1m);
  h1v=reshape(h1m,ws*ws,1);
  h2v=reshape(h2m,ws*ws,1);
  [etot(k) e2er e3er e2cm e3cm a21 a22 a31 a32]=errortot(e2v,e3v,e2m,e3m,h1v,h2v,h1m,h2m);
  k=k+1;
end
[metot K]=min(etot);
figure,plot(sig,etot);
sigma=sig(K);
sigma

%Calculating h1 using HBT Filter
h1m=hbtfilter(sigma,W);

% 0 degrees
h1m0=h1m;
h1v0=reshape(h1m0,ws*ws,1);
h2m0=imrotate(h1m0,90,'crop');
h2v0=reshape(h2m0,ws*ws,1);
[ertot(1) e2er e3er e2cm e3cm a21 a22 a31 a32]=errortot(e2v,e3v,e2m,e3m,h1v0,h2v0,h1m0,h2m0);

% 45 degrees
h1m45=imrotate(h1m,45,'crop');
h1v45=reshape(h1m45,ws*ws,1);
h2m45=imrotate(h1m45,90,'crop');
h2v45=reshape(h2m45,ws*ws,1);
[ertot(2) e2er e3er e2cm e3cm a21 a22 a31 a32]=errortot(e2v,e3v,e2m,e3m,h1v45,h2v45,h1m45,h2m45);

% 90 degrees
h1m90=imrotate(h1m,90,'crop');
h1v90=reshape(h1m90,ws*ws,1);
h2m90=imrotate(h1m90,90,'crop');
h2v90=reshape(h2m90,ws*ws,1);
[ertot(3) e2er e3er e2cm e3cm a21 a22 a31 a32]=errortot(e2v,e3v,e2m,e3m,h1v90,h2v90,h1m90,h2m90);

% 135 degrees
h1m135=imrotate(h1m,135,'crop');
h1v135=reshape(h1m135,ws*ws,1);
h2m135=imrotate(h1m135,90,'crop');
h2v135=reshape(h2m135,ws*ws,1);
[ertot(4) e2er e3er e2cm e3cm a21 a22 a31 a32]=errortot(e2v,e3v,e2m,e3m,h1v135,h2v135,h1m135,h2m135);

% Minimum error calculation and find the corresponding h1 and h2
[minertot I]=min(ertot);
switch I
    case 1 
           fh1m=h1m0;
           fh1v=h1v0;
           fh2m=h2m0;
           fh2v=h2v0;
    case 2 
           fh1m=h1m45;
           fh1v=h1v45;
           fh2m=h2m45;
           fh2v=h2v45;
    case 3 
           fh1m=h1m90;
           fh1v=h1v90;
           fh2m=h2m90;
           fh2v=h2v90;
    case 4 
           fh1m=h1m135;
           fh1v=h1v135;
           fh2m=h2m135;
           fh2v=h2v135;
end
minertot

% Gamma Calculation


 Yi=std(Ng);
 gamma=median(Yi)/0.6745;
 gamma
 c=4;
 
% Calculating Ri
 
mn=mean(Ng);
one=ones(1,ws*ws);
mnv=one'*mn;
bi=Ng-mnv+c*gamma;
Di1=sqrt(sum(bi.*bi));
Di1c=Di1';

Ni=dotprod(bi',h1v0);
nm=norm(h1v0);
Di=Di1c.*nm;
Ri(:,1)=Ni./Di;

Ni=dotprod(bi',h1v45);
nm=norm(h1v45);
Di=Di1c.*nm;
Ri(:,2)=Ni./Di;

Ni=dotprod(bi',h1v90);
nm=norm(h1v90);
Di=Di1c.*nm;
Ri(:,3)=Ni./Di;

Ni=dotprod(bi',h1v135);
nm=norm(h1v135);
Di=Di1c.*nm;
Ri(:,4)=Ni./Di;



%Largest Magnitude of similarity Maps
Ric=Ri';
Edgev=max(Ric);
Edgem=vec2mat(Edgev,row1);
Edgem=Edgem';
Edgemc=imcomplement(Edgem);
figure,imshow(Im),title('Original Image');
figure,imshow(Edgemc),title('Edge Detection');

%Applying Threshold
Threshold=0.83
EdgemT=zeros(size(Edgemc));
EdgemT=im2uint8(EdgemT);
EdgemT=EdgemT+100;
EdgemT(Edgemc>=Threshold)=255;
figure,imshow(EdgemT),title('After applying Threshold');
end


function [ertot e2er e3er e2cm e3cm a21 a22 a31 a32]=errortot(e2v,e3v,e2m,e3m,h1v,h2v,h1m,h2m)
%Calculating alpha coefficients
a21=dot(e2v,h1v)/dot(h1v,h1v);
a22=dot(e2v,h2v)/dot(h2v,h2v);
a31=dot(e3v,h1v)/dot(h1v,h1v);
a32=dot(e3v,h2v)/dot(h2v,h2v);

%Calculating approximated e2 and e3
e2cm=a21*h1m+a22*h2m;
e3cm=a31*h1m+a32*h2m;

%Calcualating errors 
e2er=norm(e2m-e2cm)/norm(e2m);
e3er=norm(e3m-e3cm)/norm(e3m);

%Calculating Total error
ertot=e2er+e3er;
end
function h1m=hbtfilter(Z,W)
for xd=-W:1:W;
    for yd=-W:1:W
        et=-Z*(xd+yd);
        h1m(xd+W+1,yd+W+1)=(1-exp(et))/(1+exp(et));
    end
end
end