function rW=denbayer1(iman,F,vr,vg,vb,scale)

%%iman:     noisy input image
%%F:        filter
%%vr,vg,vb: the levels of noises
%%scale:    the wavelet scale

%%rW:       filtered wavelet coefficients

[nf,mf]=size(F);
[n,m]=size(iman);

%%1.compute the filtered noise levels at the four positions

Fg1=zeros(nf,mf);Fr=Fg1;

%%G position 1: G R G R G R 
Fg1(1:2:nf,1:2:mf)=vg*F(1:2:nf,1:2:mf);
Fg1(2:2:nf,2:2:mf)=vg*F(2:2:nf,2:2:mf);
Fg1(1:2:nf,2:2:mf)=vr*F(1:2:nf,2:2:mf);
Fg1(2:2:nf,1:2:mf)=vb*F(2:2:nf,1:2:mf);
Vg1=sum(sum(Fg1.^2));

%%G position 2: G B G B G B 
Fg2=Fg1;
Fg2(1:2:nf,2:2:mf)=vb*F(1:2:nf,2:2:mf);
Fg2(2:2:nf,1:2:mf)=vr*F(2:2:nf,1:2:mf);
Vg2=sum(sum(Fg2.^2));

%%R position: R G R G  
Fr(1:2:nf,1:2:mf)=vr*F(1:2:nf,1:2:mf);
Fr(1:2:nf,2:2:mf)=vg*F(1:2:nf,2:2:mf);
Fr(2:2:nf,1:2:mf)=vg*F(2:2:nf,1:2:mf);
Fr(2:2:nf,2:2:mf)=vb*F(2:2:nf,2:2:mf);
Vr=sum(sum(Fr.^2));

%%B position: B G B G  
Fb=Fr;
Fb(1:2:nf,1:2:mf)=vb*F(1:2:nf,1:2:mf);
Fb(2:2:nf,2:2:mf)=vr*F(2:2:nf,2:2:mf);
Vb=sum(sum(Fb.^2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%2. process the wavelet coefficients
W=conv2(iman,F);%%convolution, size of W: (n+nf-1)x(m+mf-1)
rW=W;

%%%Note:
%   In the following denoising procedure, some boundary areas will not be processed.
%   You may extend the image boundary (e.g. using symmetrical extension) to process the full image. 

s=4*2^(scale);%%%size of the window

for i=nf:s:n+nf-1-s
    for j=mf:s:m+mf-1-s
        wd=W(i:i+s-1,j:j+s-1);%%the sxs window
        Num=s*s;
        pwd=sum(sum(wd.^2))+0.01;
        px=pwd-Num*(Vg1+Vg2+Vr+Vb)/4;
        px=max(0,px/Num);
        %%%compute the four constants
        cg1=px/(px+Vg1);
        cg2=px/(px+Vg2);
        cr=px/(px+Vr);
        cb=px/(px+Vb);
        %%%mask
        c=0.5;
        wd(1:2:s,1:2:s)=wd(1:2:s,1:2:s).*(abs(wd(1:2:s,1:2:s))>c*sqrt(Vg1));
        wd(2:2:s,2:2:s)=wd(2:2:s,2:2:s).*(abs(wd(2:2:s,2:2:s))>c*sqrt(Vg2));
        wd(1:2:s,2:2:s)=wd(1:2:s,2:2:s).*(abs(wd(1:2:s,2:2:s))>c*sqrt(Vr));
        wd(2:2:s,1:2:s)=wd(2:2:s,1:2:s).*(abs(wd(2:2:s,1:2:s))>c*sqrt(Vb));
        
        T=zeros(s,s);
        T(1:2:s,1:2:s)=cg1;
        T(2:2:s,2:2:s)=cg2;
        T(1:2:s,2:2:s)=cr;
        T(2:2:s,1:2:s)=cb;
        %%%Wiener filtering
        wd=wd.*T;
        rW(i:i+s-1,j:j+s-1)=wd;
    end
end

return;