% Treatment of boundaries.

clf

fd = fopen('barbara.raw', 'r');
X = fread(fd, [512,512], 'uchar');
X = X';
map = (0:255)'/255 * ones(1,3);
colormap(map)
imagesc(uint8(X))
axis image; set(gca,'XTick',[],'YTick',[]); title('Single stage decomposition')
title('Original image')
pause

% We will use the 9/7 filters with symmetric extension at the
% boundaries.
dwtmode('zpd')  % zpd, per, sym, sp0, sp1
%wname = 'db5'
wname = 'bior4.4'

% Compute a 2-level decomposition of the image using the 9/7 filters.
[wc,s] = wavedec2(X,2,wname);

% Extract the level 1 coefficients.
a1 = appcoef2(wc,s,wname,1);         
h1 = detcoef2('h',wc,s,1);           
v1 = detcoef2('v',wc,s,1);           
d1 = detcoef2('d',wc,s,1);           

% Extract the level 2 coefficients.
a2 = appcoef2(wc,s,wname,2);
h2 = detcoef2('h',wc,s,2);
v2 = detcoef2('v',wc,s,2);
d2 = detcoef2('d',wc,s,2);

% Display the decomposition up to level 1 only.
image(uint8([a1/2,h1*10;v1*10,d1*10]))
axis image; set(gca,'XTick',[],'YTick',[]);
title('Single stage decomposition')
pause

% Display the entire decomposition upto level 2.
st = dwtmode('status','nodisp');
if (st == 'per')
  image(uint8([[a2/4,h2*10;v2*10,d2*10],h1*10;v1*10,d1*10]));
  axis image; set(gca,'XTick',[],'YTick',[]);
  title('Two stage decomposition')
end