A=imread('Washsat.gif','gif');
N=size(A,1);
J=ceil(log2(N)); % just in case round off upwards to the nearest integer
s=sqrt(2);
e_high=80;
h=s*[0.125 0.375 0.375 0.125];
offset=-2; % actual filter coefficients are offset 
L=size(h,2);
temp1=zeros(N); % temporary container for a_{j}'s
temp2=zeros(N); % temporary container for a_{j+1}'s
temp3=zeros(N); % temporary container for d^1_{j+1}'s
temp4=zeros(N); % temporary container for d^2_{j+1}'s
temp5=zeros(N); % temporary container for the moduli
temp6=zeros(N); % temporary container for the directions
temp1=A;
scale=1;
for j=0:J-1
  for n1=1:N
    for n2=1:N
      t=0;
      for k1=1:L
        for k2=1:L
          u1=mod(n1-1+scale*(k1+offset),N)+1;
          u2=mod(n2-1+scale*(k2+offset),N)+1;
          t=t+temp1(u1,u2)*h(k1)*h(k2);
        end;
      end;
      temp2(n1,n2)=t;
      if j==0
        if n1<N
          temp3(n1,n2)=(temp1(n1+1,n2)-temp1(n1,n2))/s;
        else
          temp3(n1,n2)=(temp1(1,n2)-temp1(n1,n2))/s;
        end;
        if n2<N  
          temp4(n1,n2)=(temp1(n1,n2+1)-temp1(n1,n2))/s;
        else
          temp4(n1,n2)=(temp1(n1,1)-temp1(n1,n2))/s;
        end;
      else
        u1=mod(n1+scale/2-1,N)+1;
        u2=mod(n1-scale/2-1,N)+1;
        temp3(n1,n2)=(temp1(u1,n2)-temp1(u2,n2))/s;
        u1=mod(n2+scale/2-1,N)+1;
        u2=mod(n2-scale/2-1,N)+1;
        temp4(n1,n2)=(temp1(n1,u1)-temp1(n1,u2))/s;
      end;
      temp5(n1,n2)=sqrt(temp3(n1,n2)^2+temp4(n1,n2)^2);
      if temp4(n1,n2)==0
        temp6(n1,n2)=1; % 1 - vertical direction
      else
        p=atan(temp3(n1,n2)/temp4(n1,n2));
        if abs(p)>=3*pi/8
          temp6(n1,n2)=1; % 1 - vertical direction
        elseif p>=pi/8
          temp6(n1,n2)=2; % 2 - SW/NE direction
        elseif p<=-pi/8
          temp6(n1,n2)=4; % 4 - NW/SE direction
        else
          temp6(n1,n2)=3; % 3 - horizontal direction
        end
      end
    end;
  end;
  temp7=zeros(N); % temporary container for the modulus maxima
  for n1=2:N-1
    for n2=2:N-1
      if temp6(n1,n2)==1
        u=1;
        v=0;
      elseif temp6(n1,n2)==4
        u=-1;
        v=1;
      elseif temp6(n1,n2)==3
        u=0;
        v=1;
      else
        u=1;
        v=1;
      end
      if (temp5(n1,n2)>temp5(n1+u,n2+v))&(temp5(n1,n2)>=temp5(n1-u,n2-v))
        temp7(n1,n2)=temp5(n1,n2);
      end
%      if temp7(n1,n2)<e_high
%        temp7(n1,n2)=0;
%      end
    end
  end
  %name=strcat('horizontal_',num2str(j+1),'.bmp');
  %imwrite(uint8(temp3),name,'bmp');
  %name=strcat('vertical_',num2str(j+1),'.bmp');
  %imwrite(uint8(temp4),name,'bmp');
  t=255/max(max(temp5));
  temp5=temp5*t;
  temp7=temp7*t;
  name=strcat('moduli_',num2str(j+1),'.bmp');
  imwrite((uint8(256-temp5)),name,'bmp');
  name=strcat('moduli_maxima_',num2str(j+1),'.bmp');
  imwrite((uint8(256-temp7)),name,'bmp');
  t=255/max(max(temp6));
  temp6=temp6*t;
  name=strcat('gradient_directions_',num2str(j+1),'.bmp');
  imwrite((uint8(256-temp6)),name,'bmp');
  scale=2*scale;
  temp1=temp2;
end;

        
         
    
    
    
  