function [Lab_distance,Similarity_weight]=hvsreal2(imgo,imgc)
% Function computes a similarity measure similar to the measure as axplained in
% T. Frese, C. A. Bouman, and J. P. Allebach, "A Methodology for Designing 
% Image Similarity Metrics Based on Human Visual System Models",
% Tech. Rep. TR-ECE 97-2, Purdue University, West Lafayette, IN, (1997).

img1=double(imread(imgo,'bmp'));
img2=double(imread(imgc,'bmp'));

%    X-axis direction filter:
     Nx1=10;Sigmax1=1;Nx2=10;Sigmax2=1;Theta1=pi/2;
%    Y-axis direction filter:
     Ny1=10;Sigmay1=1;Ny2=10;Sigmay2=1;Theta2=0;

[sizex sizey dim]=size(img1);

Rimg=img1(:,:,1);
Gimg=img1(:,:,2);
Bimg=img1(:,:,3);

RCimg=img2(:,:,1);
GCimg=img2(:,:,2);
BCimg=img2(:,:,3);

for i=1:sizex
   for j=1:sizey
      x0(i,j)=0.49*Rimg(i,j)+0.31*Gimg(i,j)+0.2*Bimg(i,j);
      y0(i,j)=0.177*Rimg(i,j)+0.813*Gimg(i,j)+0.011*Bimg(i,j);
      z0(i,j)=0.01*Gimg(i,j)+0.99*Bimg(i,j);
   end
end

for i=1:sizex
   for j=1:sizey
      x0_c(i,j)=0.49*RCimg(i,j)+0.31*GCimg(i,j)+0.2*BCimg(i,j);
      y0_c(i,j)=0.177*RCimg(i,j)+0.813*GCimg(i,j)+0.011*BCimg(i,j);
      z0_c(i,j)=0.01*GCimg(i,j)+0.99*BCimg(i,j);
   end
end
1
% Perform one step decomposition.                
x1 = dwt2(x0,'db1');       
x2 = dwt2(x1,'db1');
x3 = dwt2(x2,'db1');
x4 = dwt2(x3,'db1');

y1 = dwt2(y0,'db1');       
y2 = dwt2(y1,'db1');
y3 = dwt2(y2,'db1');
y4 = dwt2(y3,'db1');

z1 = dwt2(z0,'db1');       
z2 = dwt2(z1,'db1');
z3 = dwt2(z2,'db1');
z4 = dwt2(z3,'db1');

x1_c = dwt2(x0_c,'db1');       
x2_c = dwt2(x1_c,'db1');
x3_c = dwt2(x2_c,'db1');
x4_c = dwt2(x3_c,'db1');

y1_c = dwt2(y0_c,'db1');       
y2_c = dwt2(y1_c,'db1');
y3_c = dwt2(y2_c,'db1');
y4_c = dwt2(y3_c,'db1');

z1_c = dwt2(z0_c,'db1');       
z2_c = dwt2(z1_c,'db1');
z3_c = dwt2(z2_c,'db1');
z4_c = dwt2(z3_c,'db1');
2
for i=1:sizex
   for j=1:sizey
      l0(i,j)=116*(y0(i,j).^(1/3))-16;
      a0(i,j)=500*(x0(i,j).^(1/3) - y0(i,j).^(1/3));
      b0(i,j)=200*(y0(i,j).^(1/3) - z0(i,j).^(1/3));
   end
end

for i=1:sizex/2
   for j=1:sizey/2
      l1(i,j)=116*(y1(i,j).^(1/3))-16;
      a1(i,j)=500*(x1(i,j).^(1/3) - y1(i,j).^(1/3));
      b1(i,j)=200*(y1(i,j).^(1/3) - z1(i,j).^(1/3));
   end
end

for i=1:sizex/4
   for j=1:sizey/4
      l2(i,j)=116*(y2(i,j).^(1/3))-16;
      a2(i,j)=500*(x2(i,j).^(1/3) - y2(i,j).^(1/3));
      b2(i,j)=200*(y2(i,j).^(1/3) - z2(i,j).^(1/3));
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      l3(i,j)=116*(y3(i,j).^(1/3))-16;
      a3(i,j)=500*(x3(i,j).^(1/3) - y3(i,j).^(1/3));
      b3(i,j)=200*(y3(i,j).^(1/3) - z3(i,j).^(1/3));
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      l4(i,j)=116*(y4(i,j).^(1/3))-16;
      a4(i,j)=500*(x4(i,j).^(1/3) - y4(i,j).^(1/3));
      b4(i,j)=200*(y4(i,j).^(1/3) - z4(i,j).^(1/3));
   end
end

for m=1:sizex/2
   for n=1:sizey/2
      c01y(m,n)=l0(m*2,n*2)-l1(m,n);
      c01a(m,n)=a0(m*2,n*2)-a1(m,n);
      c01b(m,n)=b0(m*2,n*2)-b1(m,n);
   end
end

for m=1:sizex/4
   for n=1:sizey/4
      c02y(m,n)=l0(m*4,n*4)-l2(m,n);
      c02a(m,n)=a0(m*4,n*4)-a2(m,n);
      c02b(m,n)=b0(m*4,n*4)-b2(m,n);
   end
end

for m=1:sizex/8
   for n=1:sizey/8
      c03y(m,n)=l0(m*8,n*8)-l3(m,n);
      c03a(m,n)=a0(m*8,n*8)-a3(m,n);
      c03b(m,n)=b0(m*8,n*8)-b3(m,n);
   end
end

for m=1:sizex/16
   for n=1:sizey/16
      c04y(m,n)=l0(m*16,n*16)-l4(m,n);
      c04a(m,n)=a0(m*16,n*16)-a4(m,n);
      c04b(m,n)=b0(m*16,n*16)-b4(m,n);
   end
end

for m=1:sizex/4
   for n=1:sizey/4
      c11y(m,n)=l1(m*2,n*2)-l2(m,n);
      c11a(m,n)=a1(m*2,n*2)-a2(m,n);
      c11b(m,n)=b1(m*2,n*2)-b2(m,n);
   end
end

for m=1:sizex/8
   for n=1:sizey/8
      c12y(m,n)=l1(m*4,n*4)-l3(m,n);
      c12a(m,n)=a1(m*4,n*4)-a3(m,n);
      c12b(m,n)=b1(m*4,n*4)-b3(m,n);
   end
end

for m=1:sizex/16
   for n=1:sizey/16
      c13y(m,n)=l1(m*8,n*8)-l4(m,n);
      c13a(m,n)=a1(m*8,n*8)-a4(m,n);
      c13b(m,n)=b1(m*8,n*8)-b4(m,n);
   end
end

for m=1:sizex/8
   for n=1:sizey/8
      c21y(m,n)=l2(m*2,n*2)-l3(m,n);
      c21a(m,n)=a2(m*2,n*2)-a3(m,n);
      c21b(m,n)=b2(m*2,n*2)-b3(m,n);
   end
end

for m=1:sizex/16
   for n=1:sizey/16
      c22y(m,n)=l2(m*4,n*4)-l4(m,n);
      c22a(m,n)=a2(m*4,n*4)-a4(m,n);
      c22b(m,n)=b2(m*4,n*4)-b4(m,n);
   end
end

for m=1:sizex/16
   for n=1:sizey/16
      c31y(m,n)=l3(m*2,n*2)-l4(m,n);
      c31a(m,n)=a3(m*2,n*2)-a4(m,n);
      c31b(m,n)=b3(m*2,n*2)-b4(m,n);
   end
end
3
T_c01y=mean2(c01y);
T_c02y=mean2(c02y);
T_c03y=mean2(c03y);
T_c04y=mean2(c04y);
T_c11y=mean2(c11y);
T_c12y=mean2(c12y);
T_c13y=mean2(c13y);
T_c21y=mean2(c21y);
T_c22y=mean2(c22y);
T_c31y=mean2(c31y);

T_c01a=mean2(c01a);
T_c02a=mean2(c02a);
T_c03a=mean2(c03a);
T_c04a=mean2(c04a);
T_c11a=mean2(c11a);
T_c12a=mean2(c12a);
T_c13a=mean2(c13a);
T_c21a=mean2(c21a);
T_c22a=mean2(c22a);
T_c31a=mean2(c31a);


T_c01b=mean2(c01b);
T_c02b=mean2(c02b);
T_c03b=mean2(c03b);
T_c04b=mean2(c04b);
T_c11b=mean2(c11b);
T_c12b=mean2(c12b);
T_c13b=mean2(c13b);
T_c21b=mean2(c21b);
T_c22b=mean2(c22b);
T_c31b=mean2(c31b);

4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:sizex
   for j=1:sizey
      l0_c(i,j)=116*(y0_c(i,j).^(1/3))-16;
      a0_c(i,j)=500*(x0_c(i,j).^(1/3) - y0_c(i,j).^(1/3));
      b0_c(i,j)=200*(y0_c(i,j).^(1/3) - z0_c(i,j).^(1/3));
   end
end

for i=1:sizex/2
   for j=1:sizey/2
      l1_c(i,j)=116*(y1_c(i,j).^(1/3))-16;
      a1_c(i,j)=500*(x1_c(i,j).^(1/3) - y1_c(i,j).^(1/3));
      b1_c(i,j)=200*(y1_c(i,j).^(1/3) - z1_c(i,j).^(1/3));
   end
end

for i=1:sizex/4
   for j=1:sizey/4
      l2_c(i,j)=116*(y2_c(i,j).^(1/3))-16;
      a2_c(i,j)=500*(x2_c(i,j).^(1/3) - y2_c(i,j).^(1/3));
      b2_c(i,j)=200*(y2_c(i,j).^(1/3) - z2_c(i,j).^(1/3));
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      l3_c(i,j)=116*(y3_c(i,j).^(1/3))-16;
      a3_c(i,j)=500*(x3_c(i,j).^(1/3) - y3_c(i,j).^(1/3));
      b3_c(i,j)=200*(y3_c(i,j).^(1/3) - z3_c(i,j).^(1/3));
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      l4_c(i,j)=116*(y4_c(i,j).^(1/3))-16;
      a4_c(i,j)=500*(x4_c(i,j).^(1/3) - y4_c(i,j).^(1/3));
      b4_c(i,j)=200*(y4_c(i,j).^(1/3) - z4_c(i,j).^(1/3));
   end
end
5
v_l0=blockimg(l0,l0_c,16,16);
v_a0=blockimg(a0,a0_c,16,16);
v_b0=blockimg(b0,b0_c,16,16);

v_l1=blockimg(l1,l1_c,8,8);
v_a1=blockimg(a1,a1_c,8,8);
v_b1=blockimg(b1,b1_c,8,8);

v_l2=blockimg(l2,l2_c,4,4);
v_a2=blockimg(a2,a2_c,4,4);
v_b2=blockimg(b2,b2_c,4,4);

v_l3=blockimg(l3,l3_c,4,4);
v_a3=blockimg(a3,a3_c,4,4);
v_b3=blockimg(b3,b3_c,4,4);
6
d_l0=mean2((l0-l0_c).^2);
d_a0=mean2((a0-a0_c).^2);
d_b0=mean2((b0-b0_c).^2);

d_l1=mean2((l1-l1_c).^2);
d_a1=mean2((a1-a1_c).^2);
d_b1=mean2((b1-b1_c).^2);

d_l2=mean2((l2-l2_c).^2);
d_a2=mean2((a2-a2_c).^2);
d_b2=mean2((b2-b2_c).^2);

d_l3=mean2((l3-l3_c).^2);
d_a3=mean2((a3-a3_c).^2);
d_b3=mean2((b3-b3_c).^2);

d_l4=mean2((l4-l4_c).^2);
d_a4=mean2((a4-a4_c).^2);
d_b4=mean2((b4-b4_c).^2);

7
for m=1:sizex/2
   for n=1:sizey/2
      c01y_c(m,n)=l0_c(m*2,n*2)-l1_c(m,n);
      c01a_c(m,n)=a0_c(m*2,n*2)-a1_c(m,n);
      c01b_c(m,n)=b0_c(m*2,n*2)-b1_c(m,n);
   end
end

for m=1:sizex/4
   for n=1:sizey/4
      c02y_c(m,n)=l0_c(m*4,n*4)-l2_c(m,n);
      c02a_c(m,n)=a0_c(m*4,n*4)-a2_c(m,n);
      c02b_c(m,n)=b0_c(m*4,n*4)-b2_c(m,n);
   end
end

for m=1:sizex/8
   for n=1:sizey/8
      c03y_c(m,n)=l0_c(m*8,n*8)-l3_c(m,n);
      c03a_c(m,n)=a0_c(m*8,n*8)-a3_c(m,n);
      c03b_c(m,n)=b0_c(m*8,n*8)-b3_c(m,n);
   end
end

for m=1:sizex/16
   for n=1:sizey/16
      c04y_c(m,n)=l0_c(m*16,n*16)-l4_c(m,n);
      c04a_c(m,n)=a0_c(m*16,n*16)-a4_c(m,n);
      c04b_c(m,n)=b0_c(m*16,n*16)-b4_c(m,n);
   end
end

for m=1:sizex/4
   for n=1:sizey/4
      c11y_c(m,n)=l1_c(m*2,n*2)-l2_c(m,n);
      c11a_c(m,n)=a1_c(m*2,n*2)-a2_c(m,n);
      c11b_c(m,n)=b1_c(m*2,n*2)-b2_c(m,n);
   end
end

for m=1:sizex/8
   for n=1:sizey/8
      c12y_c(m,n)=l1_c(m*4,n*4)-l3_c(m,n);
      c12a_c(m,n)=a1_c(m*4,n*4)-a3_c(m,n);
      c12b_c(m,n)=b1_c(m*4,n*4)-b3_c(m,n);
   end
end

for m=1:sizex/16
   for n=1:sizey/16
      c13y_c(m,n)=l1_c(m*8,n*8)-l4_c(m,n);
      c13a_c(m,n)=a1_c(m*8,n*8)-a4_c(m,n);
      c13b_c(m,n)=b1_c(m*8,n*8)-b4_c(m,n);
   end
end

for m=1:sizex/8
   for n=1:sizey/8
      c21y_c(m,n)=l2_c(m*2,n*2)-l3_c(m,n);
      c21a_c(m,n)=a2_c(m*2,n*2)-a3_c(m,n);
      c21b_c(m,n)=b2_c(m*2,n*2)-b3_c(m,n);
   end
end

for m=1:sizex/16
   for n=1:sizey/16
      c22y_c(m,n)=l2_c(m*4,n*4)-l4_c(m,n);
      c22a_c(m,n)=a2_c(m*4,n*4)-a4_c(m,n);
      c22b_c(m,n)=b2_c(m*4,n*4)-b4_c(m,n);
   end
end

for m=1:sizex/16
   for n=1:sizey/16
      c31y_c(m,n)=l3_c(m*2,n*2)-l4_c(m,n);
      c31a_c(m,n)=a3_c(m*2,n*2)-a4_c(m,n);
      c31b_c(m,n)=b3_c(m*2,n*2)-b4_c(m,n);
   end
end

8
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:sizex/2
   for j=1:sizey/2
      if c01y(i,j)>=T_c01y
         c01y(i,j)=1;
      elseif c01y(i,j)<=-T_c01y
         c01y(i,j)=-1;
      else
         c01y(i,j)=0;
      end
      if c01a(i,j)>=T_c01a
         c01a(i,j)=1;
      elseif c01a(i,j)<=-T_c01a
         c01a(i,j)=-1;
      else
         c01a(i,j)=0;
      end
      if c01b(i,j)>=T_c01b
         c01b(i,j)=1;
      elseif c01b(i,j)<=-T_c01b
         c01b(i,j)=-1;
      else
         c01b(i,j)=0;
      end
   end
end


for i=1:sizex/4
   for j=1:sizey/4
      if c02y(i,j)>=T_c02y
         c02y(i,j)=1;
      elseif c02y(i,j)<=-T_c02y
         c02y(i,j)=-1;
      else
         c02y(i,j)=0;
      end
      if c02a(i,j)>=T_c02a
         c02a(i,j)=1;
      elseif c02a(i,j)<=-T_c02a
         c02a(i,j)=-1;
      else
         c02a(i,j)=0;
      end
      if c02b(i,j)>=T_c02b
         c02b(i,j)=1;
      elseif c02b(i,j)<=-T_c02b
         c02b(i,j)=-1;
      else
         c02b(i,j)=0;
      end
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if c03y(i,j)>=T_c03y
         c03y(i,j)=1;
      elseif c03y(i,j)<=-T_c03y
         c03y(i,j)=-1;
      else
         c03y(i,j)=0;
      end
      if c03a(i,j)>=T_c03a
         c03a(i,j)=1;
      elseif c03a(i,j)<=-T_c03a
         c03a(i,j)=-1;
      else
         c03a(i,j)=0;
      end
      if c03b(i,j)>=T_c03b
         c03b(i,j)=1;
      elseif c03b(i,j)<=-T_c03b
         c03b(i,j)=-1;
      else
         c03b(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if c04y(i,j)>=T_c04y
         c04y(i,j)=1;
      elseif c04y(i,j)<=-T_c04y
         c04y(i,j)=-1;
      else
         c04y(i,j)=0;
      end
      if c04a(i,j)>=T_c04a
         c04a(i,j)=1;
      elseif c04a(i,j)<=-T_c04a
         c04a(i,j)=-1;
      else
         c04a(i,j)=0;
      end
      if c04b(i,j)>=T_c04b
         c04b(i,j)=1;
      elseif c04b(i,j)<=-T_c04b
         c04b(i,j)=-1;
      else
         c04b(i,j)=0;
      end
   end
end

for i=1:sizex/4
   for j=1:sizey/4
      if c11y(i,j)>=T_c11y
         c11y(i,j)=1;
      elseif c11y(i,j)<=-T_c11y
         c11y(i,j)=-1;
      else
         c11y(i,j)=0;
      end
      if c11a(i,j)>=T_c11a
         c11a(i,j)=1;
      elseif c11a(i,j)<=-T_c11a
         c11a(i,j)=-1;
      else
         c11a(i,j)=0;
      end
      if c11b(i,j)>=T_c11b
         c11b(i,j)=1;
      elseif c11b(i,j)<=-T_c11b
         c11b(i,j)=-1;
      else
         c11b(i,j)=0;
      end
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if c12y(i,j)>=T_c12y
         c12y(i,j)=1;
      elseif c12y(i,j)<=-T_c12y
         c12y(i,j)=-1;
      else
         c12y(i,j)=0;
      end
      if c12a(i,j)>=T_c12a
         c12a(i,j)=1;
      elseif c12a(i,j)<=-T_c12a
         c12a(i,j)=-1;
      else
         c12a(i,j)=0;
      end
      if c12b(i,j)>=T_c12b
         c12b(i,j)=1;
      elseif c12b(i,j)<=-T_c12b
         c12b(i,j)=-1;
      else
         c12b(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if c13y(i,j)>=T_c13y
         c13y(i,j)=1;
      elseif c13y(i,j)<=-T_c13y
         c13y(i,j)=-1;
      else
         c13y(i,j)=0;
      end
      if c13a(i,j)>=T_c13a
         c13a(i,j)=1;
      elseif c13a(i,j)<=-T_c13a
         c13a(i,j)=-1;
      else
         c13a(i,j)=0;
      end
      if c13b(i,j)>=T_c13b
         c13b(i,j)=1;
      elseif c13b(i,j)<=-T_c13b
         c13b(i,j)=-1;
      else
         c13b(i,j)=0;
      end
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if c21y(i,j)>=T_c21y
         c21y(i,j)=1;
      elseif c21y(i,j)<=-T_c21y
         c21y(i,j)=-1;
      else
         c21y(i,j)=0;
      end
      if c21a(i,j)>=T_c21a
         c21a(i,j)=1;
      elseif c21a(i,j)<=-T_c21a
         c21a(i,j)=-1;
      else
         c21a(i,j)=0;
      end
      if c21b(i,j)>=T_c21b
         c21b(i,j)=1;
      elseif c21b(i,j)<=-T_c21b
         c21b(i,j)=-1;
      else
         c21b(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if c22y(i,j)>=T_c22y
         c22y(i,j)=1;
      elseif c22y(i,j)<=-T_c22y
         c22y(i,j)=-1;
      else
         c22y(i,j)=0;
      end
      if c22a(i,j)>=T_c22a
         c22a(i,j)=1;
      elseif c22a(i,j)<=-T_c22a
         c22a(i,j)=-1;
      else
         c22a(i,j)=0;
      end
      if c22b(i,j)>=T_c22b
         c22b(i,j)=1;
      elseif c22b(i,j)<=-T_c22b
         c22b(i,j)=-1;
      else
         c22b(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if c31y(i,j)>=T_c31y
         c31y(i,j)=1;
      elseif c31y(i,j)<=-T_c31y
         c31y(i,j)=-1;
      else
         c31y(i,j)=0;
      end
      if c31a(i,j)>=T_c31a
         c31a(i,j)=1;
      elseif c31a(i,j)<=-T_c31a
         c31a(i,j)=-1;
      else
         c31a(i,j)=0;
      end
      if c31b(i,j)>=T_c31b
         c31b(i,j)=1;
      elseif c31b(i,j)<=-T_c31b
         c31b(i,j)=-1;
      else
         c31b(i,j)=0;
      end
   end
end
9
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:sizex/2
   for j=1:sizey/2
      if c01y_c(i,j)>=T_c01y
         c01y_c(i,j)=1;
      elseif c01y_c(i,j)<=-T_c01y
         c01y_c(i,j)=-1;
      else
         c01y_c(i,j)=0;
      end
      if c01a_c(i,j)>=T_c01a
         c01a_c(i,j)=1;
      elseif c01a_c(i,j)<=-T_c01a
         c01a_c(i,j)=-1;
      else
         c01a_c(i,j)=0;
      end
      if c01b_c(i,j)>=T_c01b
         c01b_c(i,j)=1;
      elseif c01b_c(i,j)<=-T_c01b
         c01b_c(i,j)=-1;
      else
         c01b_c(i,j)=0;
      end
   end
end


for i=1:sizex/4
   for j=1:sizey/4
      if c02y_c(i,j)>=T_c02y
         c02y_c(i,j)=1;
      elseif c02y_c(i,j)<=-T_c02y
         c02y_c(i,j)=-1;
      else
         c02y_c(i,j)=0;
      end
      if c02a_c(i,j)>=T_c02a
         c02a_c(i,j)=1;
      elseif c02a_c(i,j)<=-T_c02a
         c02a_c(i,j)=-1;
      else
         c02a_c(i,j)=0;
      end
      if c02b_c(i,j)>=T_c02b
         c02b_c(i,j)=1;
      elseif c02b_c(i,j)<=-T_c02b
         c02b_c(i,j)=-1;
      else
         c02b_c(i,j)=0;
      end
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if c03y_c(i,j)>=T_c03y
         c03y_c(i,j)=1;
      elseif c03y(i,j)<=-T_c03y
         c03y_c(i,j)=-1;
      else
         c03y_c(i,j)=0;
      end
      if c03a_c(i,j)>=T_c03a
         c03a_c(i,j)=1;
      elseif c03a_c(i,j)<=-T_c03a
         c03a_c(i,j)=-1;
      else
         c03a_c(i,j)=0;
      end
      if c03b_c(i,j)>=T_c03b
         c03b_c(i,j)=1;
      elseif c03b_c(i,j)<=-T_c03b
         c03b_c(i,j)=-1;
      else
         c03b_c(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if c04y_c(i,j)>=T_c04y
         c04y_c(i,j)=1;
      elseif c04y_c(i,j)<=-T_c04y
         c04y_c(i,j)=-1;
      else
         c04y_c(i,j)=0;
      end
      if c04a_c(i,j)>=T_c04a
         c04a_c(i,j)=1;
      elseif c04a_c(i,j)<=-T_c04a
         c04a_c(i,j)=-1;
      else
         c04a_c(i,j)=0;
      end
      if c04b_c(i,j)>=T_c04b
         c04b_c(i,j)=1;
      elseif c04b_c(i,j)<=-T_c04b
         c04b_c(i,j)=-1;
      else
         c04b_c(i,j)=0;
      end
   end
end

for i=1:sizex/4
   for j=1:sizey/4
      if c11y_c(i,j)>=T_c11y
         c11y_c(i,j)=1;
      elseif c11y_c(i,j)<=-T_c11y
         c11y_c(i,j)=-1;
      else
         c11y_c(i,j)=0;
      end
      if c11a_c(i,j)>=T_c11a
         c11a_c(i,j)=1;
      elseif c11a_c(i,j)<=-T_c11a
         c11a_c(i,j)=-1;
      else
         c11a_c(i,j)=0;
      end
      if c11b_c(i,j)>=T_c11b
         c11b_c(i,j)=1;
      elseif c11b_c(i,j)<=-T_c11b
         c11b_c(i,j)=-1;
      else
         c11b_c(i,j)=0;
      end
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if c12y_c(i,j)>=T_c12y
         c12y_c(i,j)=1;
      elseif c12y_c(i,j)<=-T_c12y
         c12y_c(i,j)=-1;
      else
         c12y_c(i,j)=0;
      end
      if c12a_c(i,j)>=T_c12a
         c12a_c(i,j)=1;
      elseif c12a_c(i,j)<=-T_c12a
         c12a_c(i,j)=-1;
      else
         c12a_c(i,j)=0;
      end
      if c12b_c(i,j)>=T_c12b
         c12b_c(i,j)=1;
      elseif c12b_c(i,j)<=-T_c12b
         c12b_c(i,j)=-1;
      else
         c12b_c(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if c13y_c(i,j)>=T_c13y
         c13y_c(i,j)=1;
      elseif c13y_c(i,j)<=-T_c13y
         c13y_c(i,j)=-1;
      else
         c13y_c(i,j)=0;
      end
      if c13a_c(i,j)>=T_c13a
         c13a_c(i,j)=1;
      elseif c13a_c(i,j)<=-T_c13a
         c13a_c(i,j)=-1;
      else
         c13a_c(i,j)=0;
      end
      if c13b_c(i,j)>=T_c13b
         c13b_c(i,j)=1;
      elseif c13b_c(i,j)<=-T_c13b
         c13b_c(i,j)=-1;
      else
         c13b_c(i,j)=0;
      end
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if c21y_c(i,j)>=T_c21y
         c21y_c(i,j)=1;
      elseif c21y_c(i,j)<=-T_c21y
         c21y_c(i,j)=-1;
      else
         c21y_c(i,j)=0;
      end
      if c21a_c(i,j)>=T_c21a
         c21a_c(i,j)=1;
      elseif c21a_c(i,j)<=-T_c21a
         c21a_c(i,j)=-1;
      else
         c21a_c(i,j)=0;
      end
      if c21b_c(i,j)>=T_c21b
         c21b_c(i,j)=1;
      elseif c21b_c(i,j)<=-T_c21b
         c21b_c(i,j)=-1;
      else
         c21b_c(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if c22y_c(i,j)>=T_c22y
         c22y_c(i,j)=1;
      elseif c22y_c(i,j)<=-T_c22y
         c22y_c(i,j)=-1;
      else
         c22y_c(i,j)=0;
      end
      if c22a_c(i,j)>=T_c22a
         c22a_c(i,j)=1;
      elseif c22a_c(i,j)<=-T_c22a
         c22a_c(i,j)=-1;
      else
         c22a_c(i,j)=0;
      end
      if c22b_c(i,j)>=T_c22b
         c22b_c(i,j)=1;
      elseif c22b_c(i,j)<=-T_c22b
         c22b_c(i,j)=-1;
      else
         c22b_c(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if c31y_c(i,j)>=T_c31y
         c31y_c(i,j)=1;
      elseif c31y_c(i,j)<=-T_c31y
         c31y_c(i,j)=-1;
      else
         c31y_c(i,j)=0;
      end
      if c31a_c(i,j)>=T_c31a
         c31a_c(i,j)=1;
      elseif c31a_c(i,j)<=-T_c31a
         c31a_c(i,j)=-1;
      else
         c31a_c(i,j)=0;
      end
      if c31b_c(i,j)>=T_c31b
         c31b_c(i,j)=1;
      elseif c31b_c(i,j)<=-T_c31b
         c31b_c(i,j)=-1;
      else
         c31b_c(i,j)=0;
      end
   end
end
10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

v_c01y=blockimg(c01y,c01y_c,16,16);
v_c01a=blockimg(c01a,c01a_c,16,16);
v_c01b=blockimg(c01b,c01b_c,16,16);

v_c02y=blockimg(c02y,c02y_c,16,16);
v_c02a=blockimg(c02a,c02a_c,16,16);
v_c02b=blockimg(c02b,c02b_c,16,16);

v_c03y=blockimg(c03y,c03y_c,16,16);
v_c03a=blockimg(c03a,c03a_c,16,16);
v_c03b=blockimg(c03b,c03b_c,16,16);

v_c04y=blockimg(c04y,c04y_c,16,16);
v_c04a=blockimg(c04a,c04a_c,16,16);
v_c04b=blockimg(c04b,c04b_c,16,16);


v_c11y=blockimg(c11y,c11y_c,8,8);
v_c11a=blockimg(c11a,c11a_c,8,8);
v_c11b=blockimg(c11b,c11b_c,8,8);

v_c12y=blockimg(c12y,c12y_c,8,8);
v_c12a=blockimg(c12a,c12a_c,8,8);
v_c12b=blockimg(c12b,c12b_c,8,8);

v_c13y=blockimg(c13y,c13y_c,8,8);
v_c13a=blockimg(c13a,c13a_c,8,8);
v_c13b=blockimg(c13b,c13b_c,8,8);


v_c21y=blockimg(c21y,c21y_c,4,4);
v_c21a=blockimg(c21a,c21a_c,4,4);
v_c21b=blockimg(c21b,c21b_c,4,4);

v_c22y=blockimg(c22y,c22y_c,4,4);
v_c22a=blockimg(c22a,c22a_c,4,4);
v_c22b=blockimg(c22b,c22b_c,4,4);

v_c31y=blockimg(c31y,c31y_c,4,4);
v_c31a=blockimg(c31a,c31a_c,4,4);
v_c31b=blockimg(c31b,c31b_c,4,4);
11
d_c01y=mean2((c01y-c01y_c).^2);
d_c01a=mean2((c01a-c01a_c).^2);
d_c01b=mean2((c01b-c01b_c).^2);

d_c02y=mean2((c02y-c02y_c).^2);
d_c02a=mean2((c02a-c02a_c).^2);
d_c02b=mean2((c02b-c02b_c).^2);

d_c03y=mean2((c03y-c03y_c).^2);
d_c03a=mean2((c03a-c03a_c).^2);
d_c03b=mean2((c03b-c03b_c).^2);

d_c04y=mean2((c04y-c04y_c).^2);
d_c04a=mean2((c04a-c04a_c).^2);
d_c04b=mean2((c04b-c04b_c).^2);


d_c11y=mean2((c11y-c11y_c).^2);
d_c11a=mean2((c11a-c11a_c).^2);
d_c11b=mean2((c11b-c11b_c).^2);

d_c12y=mean2((c12y-c12y_c).^2);
d_c12a=mean2((c12a-c12a_c).^2);
d_c12b=mean2((c12b-c12b_c).^2);

d_c13y=mean2((c13y-c13y_c).^2);
d_c13a=mean2((c13a-c13a_c).^2);
d_c13b=mean2((c13b-c13b_c).^2);

d_c21y=mean2((c21y-c21y_c).^2);
d_c21a=mean2((c21a-c21a_c).^2);
d_c21b=mean2((c21b-c21b_c).^2);

d_c22y=mean2((c22y-c22y_c).^2);
d_c22a=mean2((c22a-c22a_c).^2);
d_c22b=mean2((c22b-c22b_c).^2);

d_c31y=mean2((c31y-c31y_c).^2);
d_c31a=mean2((c31a-c31a_c).^2);
d_c31b=mean2((c31b-c31b_c).^2);

12

Lab_distance= mean2((l0-l0_c).^2 +(a0-a0_c).^2 +(b0-b0_c).^2); 
Lab_distance=sqrt(Lab_distance)


filterx=d2dgauss(Nx1,Sigmax1,Nx2,Sigmax2,Theta1);
l0x= conv2(l0,filterx,'same');
a0x= conv2(a0,filterx,'same');
b0x= conv2(b0,filterx,'same');

l1x= conv2(l1,filterx,'same');
a1x= conv2(a1,filterx,'same');
b1x= conv2(b1,filterx,'same');

l2x= conv2(l2,filterx,'same');
a2x= conv2(a2,filterx,'same');
b2x= conv2(b2,filterx,'same');

l3x= conv2(l3,filterx,'same');
a3x= conv2(a3,filterx,'same');
b3x= conv2(b3,filterx,'same');

l4x= conv2(l4,filterx,'same');
a4x= conv2(a4,filterx,'same');
b4x= conv2(b4,filterx,'same');

13

filtery=d2dgauss(Ny1,Sigmay1,Ny2,Sigmay2,Theta2);

l0y= conv2(l0,filtery,'same');
a0y= conv2(a0,filtery,'same');
b0y= conv2(b0,filtery,'same');

l1y= conv2(l1,filtery,'same');
a1y= conv2(a1,filtery,'same');
b1y= conv2(b1,filtery,'same');

l2y= conv2(l2,filtery,'same');
a2y= conv2(a2,filtery,'same');
b2y= conv2(b2,filtery,'same');

l3y= conv2(l3,filtery,'same');
a3y= conv2(a3,filtery,'same');
b3y= conv2(b3,filtery,'same');

l4y= conv2(l4,filtery,'same');
a4y= conv2(a4,filtery,'same');
b4y= conv2(b4,filtery,'same');
14

% Norm of the gradient (Combining the X and Y directional derivatives)
vl0=atan2(l0y,l0x);
va0=atan2(a0y,a0x);
vb0=atan2(b0y,b0x);

vl1=atan2(l1y,l1x);
va1=atan2(a1y,a1x);
vb1=atan2(b1y,b1x);

vl2=atan2(l2y,l2x);
va2=atan2(a2y,a2x);
vb2=atan2(b2y,b2x);

vl3=atan2(l3y,l3x);
va3=atan2(a3y,a3x);
vb3=atan2(b3y,b3x);

vl4=atan2(l4y,l4x);
va4=atan2(a4y,a4x);
vb4=atan2(b4y,b4x);


sl0=sqrt(l0x.*l0x+l0y.*l0y);
sa0=sqrt(a0x.*a0x+a0y.*a0y);
sb0=sqrt(b0x.*b0x+b0y.*b0y);

sl1=sqrt(l1x.*l1x+l1y.*l1y);
sa1=sqrt(a1x.*a1x+a1y.*a1y);
sb1=sqrt(b1x.*b1x+b1y.*b1y);

sl2=sqrt(l2x.*l2x+l2y.*l2y);
sa2=sqrt(a2x.*a2x+a2y.*a2y);
sb2=sqrt(b2x.*b2x+b2y.*b2y);

sl3=sqrt(l3x.*l3x+l3y.*l3y);
sa3=sqrt(a3x.*a3x+a3y.*a3y);
sb3=sqrt(b3x.*b3x+b3y.*b3y);

sl4=sqrt(l4x.*l4x+l4y.*l4y);
sa4=sqrt(a4x.*a4x+a4y.*a4y);
sb4=sqrt(b4x.*b4x+b4y.*b4y);
15
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

l0x_c= conv2(l0_c,filterx,'same');
a0x_c= conv2(a0_c,filterx,'same');
b0x_c= conv2(b0_c,filterx,'same');

l1x_c= conv2(l1_c,filterx,'same');
a1x_c= conv2(a1_c,filterx,'same');
b1x_c= conv2(b1_c,filterx,'same');

l2x_c= conv2(l2_c,filterx,'same');
a2x_c= conv2(a2_c,filterx,'same');
b2x_c= conv2(b2_c,filterx,'same');

l3x_c= conv2(l3_c,filterx,'same');
a3x_c= conv2(a3_c,filterx,'same');
b3x_c= conv2(b3_c,filterx,'same');

l4x_c= conv2(l4_c,filterx,'same');
a4x_c= conv2(a4_c,filterx,'same');
b4x_c= conv2(b4_c,filterx,'same');


l0y_c= conv2(l0_c,filtery,'same');
a0y_c= conv2(a0_c,filtery,'same');
b0y_c= conv2(b0_c,filtery,'same');

l1y_c= conv2(l1_c,filtery,'same');
a1y_c= conv2(a1_c,filtery,'same');
b1y_c= conv2(b1_c,filtery,'same');

l2y_c= conv2(l2_c,filtery,'same');
a2y_c= conv2(a2_c,filtery,'same');
b2y_c= conv2(b2_c,filtery,'same');

l3y_c= conv2(l3_c,filtery,'same');
a3y_c= conv2(a3_c,filtery,'same');
b3y_c= conv2(b3_c,filtery,'same');

l4y_c= conv2(l4_c,filtery,'same');
a4y_c= conv2(a4_c,filtery,'same');
b4y_c= conv2(b4_c,filtery,'same');


% Norm of the gradient (Combining the X and Y directional derivatives)
vl0_c=atan2(l0y_c,l0x_c);
va0_c=atan2(a0y_c,a0x_c);
vb0_c=atan2(b0y_c,b0x_c);

vl1_c=atan2(l1y_c,l1x_c);
va1_c=atan2(a1y_c,a1x_c);
vb1_c=atan2(b1y_c,b1x_c);

vl2_c=atan2(l2y_c,l2x_c);
va2_c=atan2(a2y_c,a2x_c);
vb2_c=atan2(b2y_c,b2x_c);

vl3_c=atan2(l3y_c,l3x_c);
va3_c=atan2(a3y_c,a3x_c);
vb3_c=atan2(b3y_c,b3x_c);

vl4_c=atan2(l4y_c,l4x_c);
va4_c=atan2(a4y_c,a4x_c);
vb4_c=atan2(b4y_c,b4x_c);


sl0_c=sqrt(l0x_c.*l0x_c+l0y_c.*l0y_c);
sa0_c=sqrt(a0x_c.*a0x_c+a0y_c.*a0y_c);
sb0_c=sqrt(b0x_c.*b0x_c+b0y_c.*b0y_c);

sl1_c=sqrt(l1x_c.*l1x_c+l1y_c.*l1y_c);
sa1_c=sqrt(a1x_c.*a1x_c+a1y_c.*a1y_c);
sb1_c=sqrt(b1x_c.*b1x_c+b1y_c.*b1y_c);

sl2_c=sqrt(l2x_c.*l2x_c+l2y_c.*l2y_c);
sa2_c=sqrt(a2x_c.*a2x_c+a2y_c.*a2y_c);
sb2_c=sqrt(b2x_c.*b2x_c+b2y_c.*b2y_c);

sl3_c=sqrt(l3x_c.*l3x_c+l3y_c.*l3y_c);
sa3_c=sqrt(a3x_c.*a3x_c+a3y_c.*a3y_c);
sb3_c=sqrt(b3x_c.*b3x_c+b3y_c.*b3y_c);

sl4_c=sqrt(l4x_c.*l4x_c+l4y_c.*l4y_c);
sa4_c=sqrt(a4x_c.*a4x_c+a4y_c.*a4y_c);
sb4_c=sqrt(b4x_c.*b4x_c+b4y_c.*b4y_c);

16
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T_sl0=mean2(sl0);
T_sa0=mean2(sa0);
T_sb0=mean2(sb0);

T_sl1=mean2(sl1);
T_sa1=mean2(sa1);
T_sb1=mean2(sb1);

T_sl2=mean2(sl2);
T_sa2=mean2(sa2);
T_sb2=mean2(sb2);

T_sl3=mean2(sl3);
T_sa3=mean2(sa3);
T_sb3=mean2(sb3);

T_sl4=mean2(sl4);
T_sa4=mean2(sa4);
T_sb4=mean2(sb4);


17
for i=1:sizex
   for j=1:sizey
      if sl0(i,j)>=T_sl0
         sl0(i,j)=1;
      else
         sl0(i,j)=0;
      end
      if sa0(i,j)>=T_sa0
         sa0(i,j)=1;
      else
         sa0(i,j)=0;
      end
      if sb0(i,j)>=T_sb0
         sb0(i,j)=1;
      else
         sb0(i,j)=0;
      end
   end
end

for i=1:sizex/2
   for j=1:sizey/2
      if sl1(i,j)>=T_sl1
         sl1(i,j)=1;
      else
         sl1(i,j)=0;
      end
      if sa1(i,j)>=T_sa1
         sa1(i,j)=1;
      else
         sa1(i,j)=0;
      end
      if sb1(i,j)>=T_sb1
         sb1(i,j)=1;
      else
         sb1(i,j)=0;
      end
   end
end

for i=1:sizex/4
   for j=1:sizey/4
      if sl2(i,j)>=T_sl2
         sl2(i,j)=1;
      else
         sl2(i,j)=0;
      end
      if sa2(i,j)>=T_sa2
         sa2(i,j)=1;
      else
         sa2(i,j)=0;
      end
      if sb2(i,j)>=T_sb2
         sb2(i,j)=1;
      else
         sb2(i,j)=0;
      end
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if sl3(i,j)>=T_sl3
         sl3(i,j)=1;
      else
         sl3(i,j)=0;
      end
      if sa3(i,j)>=T_sa3
         sa3(i,j)=1;
      else
         sa3(i,j)=0;
      end
      if sb3(i,j)>=T_sb3
         sb3(i,j)=1;
      else
         sb3(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if sl4(i,j)>=T_sl4
         sl4(i,j)=1;
      else
         sl4(i,j)=0;
      end
      if sa4(i,j)>=T_sa4
         sa4(i,j)=1;
      else
         sa4(i,j)=0;
      end
      if sb4(i,j)>=T_sb4
         sb4(i,j)=1;
      else
         sb4(i,j)=0;
      end
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:sizex
   for j=1:sizey
      if sl0_c(i,j)>=T_sl0
         sl0_c(i,j)=1;
      else
         sl0_c(i,j)=0;
      end
      if sa0_c(i,j)>=T_sa0
         sa0_c(i,j)=1;
      else
         sa0_c(i,j)=0;
      end
      if sb0_c(i,j)>=T_sb0
         sb0_c(i,j)=1;
      else
         sb0_c(i,j)=0;
      end
   end
end

for i=1:sizex/2
   for j=1:sizey/2
      if sl1_c(i,j)>=T_sl1
         sl1_c(i,j)=1;
      else
         sl1_c(i,j)=0;
      end
      if sa1_c(i,j)>=T_sa1
         sa1_c(i,j)=1;
      else
         sa1_c(i,j)=0;
      end
      if sb1_c(i,j)>=T_sb1
         sb1_c(i,j)=1;
      else
         sb1_c(i,j)=0;
      end
   end
end

for i=1:sizex/4
   for j=1:sizey/4
      if sl2_c(i,j)>=T_sl2
         sl2_c(i,j)=1;
      else
         sl2_c(i,j)=0;
      end
      if sa2_c(i,j)>=T_sa2
         sa2_c(i,j)=1;
      else
         sa2_c(i,j)=0;
      end
      if sb2_c(i,j)>=T_sb2
         sb2_c(i,j)=1;
      else
         sb2_c(i,j)=0;
      end
   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if sl3_c(i,j)>=T_sl3
         sl3_c(i,j)=1;
      else
         sl3_c(i,j)=0;
      end
      if sa3_c(i,j)>=T_sa3
         sa3_c(i,j)=1;
      else
         sa3_c(i,j)=0;
      end
      if sb3_c(i,j)>=T_sb3
         sb3_c(i,j)=1;
      else
         sb3_c(i,j)=0;
      end
   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if sl4_c(i,j)>=T_sl4
         sl4_c(i,j)=1;
      else
         sl4_c(i,j)=0;
      end
      if sa4_c(i,j)>=T_sa4
         sa4_c(i,j)=1;
      else
         sa4_c(i,j)=0;
      end
      if sb4_c(i,j)>=T_sb4
         sb4_c(i,j)=1;
      else
         sb4_c(i,j)=0;
      end
   end
end

18

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:sizex
   for j=1:sizey
      if sl0(i,j)==1 & sl0_c(i,j)==1
         dvl0(i,j)=rem( (vl0(i,j)- vl0_c(i,j)),2*pi );
      elseif sl0(i,j)==1 & sl0_c(i,j)==0
         dvl0(i,j)=pi;
      elseif sl0(i,j)==0
         dvl0(i,j)=0;
      end
      
      if sa0(i,j)==1 & sa0_c(i,j)==1
         dva0(i,j)=rem( (va0(i,j)- va0_c(i,j)),2*pi );
      elseif sa0(i,j)==1 & sa0_c(i,j)==0
         dva0(i,j)=pi;
      elseif sa0(i,j)==0
         dva0(i,j)=0;
      end
      
      if sb0(i,j)==1 & sb0_c(i,j)==1
         dvb0(i,j)=rem( (vb0(i,j)- vb0_c(i,j)),2*pi );
      elseif sb0(i,j)==1 & sb0_c(i,j)==0
         dvb0(i,j)=pi;
      elseif sb0(i,j)==0
         dvb0(i,j)=0;
      end

   end
end

for i=1:sizex/2
   for j=1:sizey/2
      if sl1(i,j)==1 & sl1_c(i,j)==1
         dvl1(i,j)=rem( (vl1(i,j)- vl1_c(i,j)),2*pi );
      elseif sl1(i,j)==1 & sl1_c(i,j)==0
         dvl1(i,j)=pi;
      elseif sl1(i,j)==0
         dvl1(i,j)=0;
      end
      
      if sa1(i,j)==1 & sa1_c(i,j)==1
         dva1(i,j)=rem( (va1(i,j)- va1_c(i,j)),2*pi );
      elseif sa1(i,j)==1 & sa1_c(i,j)==0
         dva1(i,j)=pi;
      elseif sa1(i,j)==0
         dva1(i,j)=0;
      end
      
      if sb1(i,j)==1 & sb1_c(i,j)==1
         dvb1(i,j)=rem( (vb1(i,j)- vb1_c(i,j)),2*pi );
      elseif sb1(i,j)==1 & sb1_c(i,j)==0
         dvb1(i,j)=pi;
      elseif sb1(i,j)==0
         dvb1(i,j)=0;
      end

   end
end

for i=1:sizex/4
   for j=1:sizey/4
      if sl2(i,j)==1 & sl2_c(i,j)==1
         dvl2(i,j)=rem( (vl2(i,j)- vl2_c(i,j)),2*pi );
      elseif sl2(i,j)==1 & sl2_c(i,j)==0
         dvl2(i,j)=pi;
      elseif sl2(i,j)==0
         dvl2(i,j)=0;
      end
      
      if sa2(i,j)==1 & sa2_c(i,j)==1
         dva2(i,j)=rem( (va2(i,j)- va2_c(i,j)),2*pi );
      elseif sa2(i,j)==1 & sa2_c(i,j)==0
         dva2(i,j)=pi;
      elseif sa2(i,j)==0
         dva2(i,j)=0;
      end
      
      if sb2(i,j)==1 & sb2_c(i,j)==1
         dvb2(i,j)=rem( (vb2(i,j)- vb2_c(i,j)),2*pi );
      elseif sb2(i,j)==1 & sb2_c(i,j)==0
         dvb2(i,j)=pi;
      elseif sb2(i,j)==0
         dvb2(i,j)=0;
      end

   end
end

for i=1:sizex/8
   for j=1:sizey/8
      if sl3(i,j)==1 & sl3_c(i,j)==1
         dvl3(i,j)=rem( (vl3(i,j)- vl3_c(i,j)),2*pi );
      elseif sl3(i,j)==1 & sl3_c(i,j)==0
         dvl3(i,j)=pi;
      elseif sl3(i,j)==0
         dvl3(i,j)=0;
      end
      
      if sa3(i,j)==1 & sa3_c(i,j)==1
         dva3(i,j)=rem( (va3(i,j)- va3_c(i,j)),2*pi );
      elseif sa3(i,j)==1 & sa3_c(i,j)==0
         dva3(i,j)=pi;
      elseif sa3(i,j)==0
         dva3(i,j)=0;
      end
      
      if sb3(i,j)==1 & sb3_c(i,j)==1
         dvb3(i,j)=rem( (vb3(i,j)- vb3_c(i,j)),2*pi );
      elseif sb3(i,j)==1 & sb3_c(i,j)==0
         dvb3(i,j)=pi;
      elseif sb3(i,j)==0
         dvb3(i,j)=0;
      end

   end
end

for i=1:sizex/16
   for j=1:sizey/16
      if sl4(i,j)==1 & sl4_c(i,j)==1
         dvl4(i,j)=rem( (vl4(i,j)- vl4_c(i,j)),2*pi );
      elseif sl4(i,j)==1 & sl4_c(i,j)==0
         dvl4(i,j)=pi;
      elseif sl4(i,j)==0
         dvl4(i,j)=0;
      end
      
      if sa4(i,j)==1 & sa4_c(i,j)==1
         dva4(i,j)=rem( (va4(i,j)- va4_c(i,j)),2*pi );
      elseif sa4(i,j)==1 & sa4_c(i,j)==0
         dva4(i,j)=pi;
      elseif sa4(i,j)==0
         dva4(i,j)=0;
      end
      
      if sb4(i,j)==1 & sb4_c(i,j)==1
         dvb4(i,j)=rem( (vb4(i,j)- vb4_c(i,j)),2*pi );
      elseif sb4(i,j)==1 & sb4_c(i,j)==0
         dvb4(i,j)=pi;
      elseif sb4(i,j)==0
         dvb4(i,j)=0;
      end

   end
end
19
thtl0=mean2(dvl0.^2);
thta0=mean2(dva0.^2);
thtb0=mean2(dvb0.^2);

thtl1=mean2(dvl1.^2);
thta1=mean2(dva1.^2);
thtb1=mean2(dvb1.^2);

thtl2=mean2(dvl2.^2);
thta2=mean2(dva2.^2);
thtb2=mean2(dvb2.^2);

thtl3=mean2(dvl3.^2);
thta3=mean2(dva3.^2);
thtb3=mean2(dvb3.^2);

thtl4=mean2(dvl4.^2);
thta4=mean2(dva4.^2);
thtb4=mean2(dvb4.^2);
      20
thtlab=0;         
   disc=0;         
   varc=0;
   varlab=0;
   dislab=0;
   Similarity_weight=0.177*d_b3 +0.599*d_l3+ 0.962*d_a4 + 1*d_b4+ 0.289*v_c01y + 0.223*v_c01a + 0.230*v_c11y + 0.256*v_c21a +  0.148*thtl2 + 0.286*thtl3 + 0.194*thtb3 + 0.972*thta4 + 0.422*thtb4