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 denoisinghh1 = daubcqf(16);hh2 = daubcqf(8);hh3 = daubcqf(10);hh4 = daubcqf(12);hh5 = daubcqf(14);%wavelet trasnformswxn1 = mdwt(xn,hh1);pp= wxn1;%generate pilot estimate using usual wavelet denoising%and then estimate radius of Besov ballstd = 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 routinefor z = 1:ni%first project onto Besov ball in hh1 wavelet domainfor 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 domainpp = mdwt(midwt(p,hh1),hh2);%project in hh2 domainfor 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 domainpp = mdwt(midwt(p,hh2),hh3);%project in hh3 domainfor 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 domainfor 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 estimatey = midwt(p,hh5);