besovdenoise.m

You need to download the general purpose toolbox and the signal toolbox. You need to unzip these

function y=besovdenoise(xn,nr,ni)
%function y=besovdenoise(xn,nr,ni)
%
%Multiple Wavelet Basis Denoising using Besov Projections
%Requires Rice Wavelet Toolbox (www.dsp.rice.edu/software)
%
% xn: noisy image
%     works best when pixel values of original image are in [0,1]
% nr: factor used to estimate Besov norm
%     recommended value: 2.0 - 3.0
% ni: number of POCS iteration
%     recommented value: 10
% y: denoised image
%
%
% functions called:
% b_11i_norm2.m, l1_pr2.m, l1_pr.m
%
% Last modified October 30, 2003
%

%wavelet bases used for denoising
hh1 = daubcqf(16);
hh2 = daubcqf(8);
hh3 = daubcqf(10);
hh4 = daubcqf(12);
hh5 = daubcqf(14);

%wavelet trasnforms
wxn1 = mdwt(xn,hh1);

pp= wxn1;

%generate pilot estimate using usual wavelet denoising
%and then estimate radius of Besov balls

td = denoise(xn,hh1,0,[0 2.23 0 0 0 0]);
r1 = nr*b_11i_norm2(mdwt(td,hh1));
td = denoise(xn,hh2,0,[0 2.23 0 0 0 0]);
r2 = nr*b_11i_norm2(mdwt(td,hh2));
td = denoise(xn,hh3,0,[0 2.23 0 0 0 0]);
r3 = nr*b_11i_norm2(mdwt(td,hh3));
td = denoise(xn,hh4,0,[0 2.23 0 0 0 0]);
r4 = nr*b_11i_norm2(mdwt(td,hh4));
td = denoise(xn,hh5,0,[0 2.23 0 0 0 0]);
r5 = nr*b_11i_norm2(mdwt(td,hh5));

clear td;

[m,n]=size(xn);
scale = ones(size(xn));

level = log2(n);

%main POCS routine
for z = 1:ni

%first project onto Besov ball in hh1 wavelet domain
for k=1:level
 ei = 2^k;
 si = ei/2 + 1;
 bhh = pp(si:ei,si:ei); ahh = scale(si:ei,si:ei);
 bhl = pp(1:si-1,si:ei); ahl = scale(1:si-1,si:ei);
 blh = pp(si:ei,1:si-1); alh = scale(si:ei,1:si-1);

 temp = sum(sum(abs(bhh).*ahh+abs(bhl).*ahl+abs(blh).*alh))+pp(1,1); 
 if (temp > r1)
  [phh,phl,plh]=l1_pr2(bhh,bhl,blh,ahh,ahl,alh,r1-abs(pp(1,1)));
  p(si:ei,si:ei)=phh;
  p(1:si-1,si:ei)=phl;
  p(si:ei,1:si-1)=plh;
 else
  p(si:ei,si:ei)=bhh;
  p(1:si-1,si:ei)=bhl;
  p(si:ei,1:si-1)=blh;
 end;
end;
p(1) = pp(1);

%wavelet transform of hh1 projected image in hh2 domain
pp = mdwt(midwt(p,hh1),hh2);

%project in hh2 domain
for k=1:level
 ei = 2^k;
 si = ei/2 + 1;
 bhh = pp(si:ei,si:ei); ahh = scale(si:ei,si:ei);
 bhl = pp(1:si-1,si:ei); ahl = scale(1:si-1,si:ei);
 blh = pp(si:ei,1:si-1); alh = scale(si:ei,1:si-1);
 
 temp = sum(sum(abs(bhh).*ahh+abs(bhl).*ahl+abs(blh).*alh))+pp(1,1);
 if (temp > r2)
  [phh,phl,plh]=l1_pr2(bhh,bhl,blh,ahh,ahl,alh,r2-abs(pp(1,1)));
  p(si:ei,si:ei)=phh;
  p(1:si-1,si:ei)=phl;
  p(si:ei,1:si-1)=plh;
 else
  p(si:ei,si:ei)=bhh;
  p(1:si-1,si:ei)=bhl;
  p(si:ei,1:si-1)=blh;
 end;
end;
p(1) = pp(1);

%wavelet transform of hh2 projected image in hh3 domain
pp = mdwt(midwt(p,hh2),hh3);

%project in hh3 domain
for k=1:level
 ei = 2^k;
 si = ei/2 + 1;
 bhh = pp(si:ei,si:ei); ahh = scale(si:ei,si:ei);
 bhl = pp(1:si-1,si:ei); ahl = scale(1:si-1,si:ei);
 blh = pp(si:ei,1:si-1); alh = scale(si:ei,1:si-1);
 
 temp = sum(sum(abs(bhh).*ahh+abs(bhl).*ahl+abs(blh).*alh))+pp(1,1);
 if (temp > r3)
  [phh,phl,plh]=l1_pr2(bhh,bhl,blh,ahh,ahl,alh,r3-abs(pp(1,1)));
  p(si:ei,si:ei)=phh;
  p(1:si-1,si:ei)=phl;
  p(si:ei,1:si-1)=plh;
 else
  p(si:ei,si:ei)=bhh;
  p(1:si-1,si:ei)=bhl;
  p(si:ei,1:si-1)=blh;
 end;
end;
p(1) = pp(1);
 

pp = mdwt(midwt(p,hh3),hh4);

%now project in hh4 domain 
for k=1:level
 ei = 2^k;
 si = ei/2 + 1;
 bhh = pp(si:ei,si:ei); ahh = scale(si:ei,si:ei);
 bhl = pp(1:si-1,si:ei); ahl = scale(1:si-1,si:ei);
 blh = pp(si:ei,1:si-1); alh = scale(si:ei,1:si-1);
 
 temp = sum(sum(abs(bhh).*ahh+abs(bhl).*ahl+abs(blh).*alh))+pp(1,1);
 if (temp > r4)
  [phh,phl,plh]=l1_pr2(bhh,bhl,blh,ahh,ahl,alh,r4-abs(pp(1,1)));
  p(si:ei,si:ei)=phh;
  p(1:si-1,si:ei)=phl;
  p(si:ei,1:si-1)=plh;
 else
  p(si:ei,si:ei)=bhh;
  p(1:si-1,si:ei)=bhl;
  p(si:ei,1:si-1)=blh;
 end;
end;
p(1) = pp(1);
 
pp = mdwt(midwt(p,hh4),hh5);

%projection in hh5 domain
for k=1:level
 ei = 2^k;
 si = ei/2 + 1;
 bhh = pp(si:ei,si:ei); ahh = scale(si:ei,si:ei);
 bhl = pp(1:si-1,si:ei); ahl = scale(1:si-1,si:ei);
 blh = pp(si:ei,1:si-1); alh = scale(si:ei,1:si-1);
 
 temp = sum(sum(abs(bhh).*ahh+abs(bhl).*ahl+abs(blh).*alh))+pp(1,1);
 if (temp > r4)
  [phh,phl,plh]=l1_pr2(bhh,bhl,blh,ahh,ahl,alh,r4-abs(pp(1,1)));
  p(si:ei,si:ei)=phh;
  p(1:si-1,si:ei)=phl;
  p(si:ei,1:si-1)=plh;
 else
  p(si:ei,si:ei)=bhh;
  p(1:si-1,si:ei)=bhl;
  p(si:ei,1:si-1)=blh;
 end;
end;
p(1) = pp(1);

%bring it back to hh1 domain 
pp = mdwt(midwt(p,hh5),hh1);

end; %of z (iteration)

%invert wavelet transform to get final estimate
y = midwt(p,hh5);