%%Read 512x512 image in
clear all
colormap(gray(256))
pic=imread('bsb.tif');
pic=double(pic);

colormap(gray(256))
%Break pic into 256 (32x32) blocks
intermimage=zeros(16,16);
finalimage=zeros(512,512);
blocknum=0;
vcoeffmatrix=zeros(256,512);
xc=0;
yc=0;
varmatrix=zeros(16,16);
for I = 1 : 32: 481,
  xc=xc+1;
  for J= 1 : 32: 481,
    yc=yc+1;
    blocknum=blocknum+1;
    block= pic(I:I+31,J:J+31);
    h=daubcqf(4,'min');
    L=1;
    [w,L]=mdwt(block,h,L); 

    %%w is wavelet transform of 32x32 block
    LH=w(17:32,1:16);        %stores vertical edge detector coeffs in matrix LH
    HL=w(1:16,17:32);        %stores horizontal edge detector coeffs in matrix HL


    v32=LH(1,:);               %store wavelet coeffs for LH in first half of vector v
    for k = 2 : 16,
      v32=[v32,LH(k,:)];
    end
    for k= 1 :16,            %store wavelet coeffs for HL in second half of vector v
      v32=[v32,HL(k,:)];
    end

    vcoeffmatrix(blocknum,:)=v32;
      
    variance=var(v32);         %variance of LH and HL wavelet coeffs
    varmatrix(xc,yc)=variance;
    
    classify512
    %%classifies 32x32 block in intermediate image

  end
  yc=0;
end

%%Intermimage is the 16x16 intermediate image, containing 32x32 blocks, some of which may
%%be undetermined

blocknum=0;
for I=1:16,
  for J=1:16,
    blocknum=blocknum+1;   
    makefinal
    %%classifies 32x32 block in final image, also classifies undetermined
    %%32x32 blocks to produce the final image
  end
end


figure(1)
imagesc(pic)
title('Original Image')
figure(2)
imagesc(finalimage)
title('Segmented Image')
