denoise_band.m

MATLAB Denoising software for grayscale images

function newBand = denoise_band(pyr,pyrN,pind,nband,opts,prnt);

% function newBand = denoise_band(pyr,pyrN,pind,nband,opts,parent);
%
% Implementation of the parameter estimation and iterated wiener filter
%
% where
%   'pyr'       is the wavelet decompostion of the noisy band
%   'pyrN'      WP of the noise to estimate covariance
%   'pind'      indices of the subbands within 'pyr'
%   'nband'     number of the band to denoise
%   'opts'      set of options, see 'newOptions'
%   'prnt'      1 if parent coefficient should be included, 0 otherwise
%
%   'newBand'   denoised version of the subband 
%
% Peter-Vincent Gehler 31.May 2005

verbose = 0;

% ... get the current bands ...
band1 = pyrBand(pyr,pind,nband);
band1n = pyrBand(pyrN,pind,nband);

% ... and its parent band if needed...
if prnt
    band2 = pyrBand(pyr,pind,nband+opts.nOrientations);
    band2n = pyrBand(pyrN,pind,nband+opts.nOrientations) ;
    if sum(size(band2)<size(band1))
        band2 = real(expand(band2,2));
        band2n = real(expand(band2n,2));
    end
    band2n = band2n;
else
    band2 = []; band2n = [];
end

if strcmp(opts.pyrtype,'fullsf') | strcmp(opts.pyrtype,'sf')
    [Nsy,Nsx] = size(band1n);
    band1n = band1n * sqrt(((Nsy-2)*(Nsx-2))/(Nsy*Nsx));
    band2n = band2n * sqrt(((Nsy-2)*(Nsx-2))/(Nsy*Nsx));
end

% ... the pyramid representations are not longer needed, free the mem ...
clear pyr pyrN

% The subband neighbourhood adaption ...
if opts.adaptneighbours
    nblv = size(band1,1)-opts.bigblock(1)+1;
    nblh = size(band1,2)-opts.bigblock(2)+1;
    Lx = (opts.bigblock(1)-1)/2;
    Ly = (opts.bigblock(2)-1)/2;
    
    % ... get all kxk neighbourhoods from the band, from these
    % coefficients the neighbourhood is chosen ... 
    [Y,center] = extractNeighbourhood(band1,opts.bigblock,band2,prnt);
    nCoeff = size(Y,1);
    
    % ... this function computes the cross cumulant, in 'index' the
    % indices are stored to access the neighbourhood within the
    % larger kxk neighbourhood, D is the number of coefficients used ...
    [index,D] = adaptiveNeighbourhood(Y(:,1:end-prnt)',center,opts,prnt);
    
    center = 1;
    
    % ... and keep only desired neighbours ...
    Y = Y(:,index);
    
    % plot something 
    if verbose > 1
        fprintf('neighbourhood size %d\n',D-prnt);
    end
    if verbose > 2
        TT = zeros(opts.bigblock);
        TT(index(1:end-prnt)) = 1;
        myimagesc(TT,3);
        title('neighbourhoud used');
        drawnow;
        fprintf('will use a neighbourhood of %d coeffs',D-prnt);
        if prnt,fprintf(' plus parent\n');
        else,fprintf('\n');
        end
    end
else
    % ... will use a 'block' surrounding the center pixel ...
    D = prod(opts.block) + prnt;
    nblv = size(band1,1)-opts.block(1)+1;
    nblh = size(band1,2)-opts.block(2)+1;
    nCoeff = (size(band1,1)-opts.block(1)+1) * (size(band1,2)-opts.block(2)+1);
    Lx = (opts.block(1)-1)/2;
    Ly = (opts.block(2)-1)/2;
    [Y,center] = extractNeighbourhood(band1,opts.block,band2,prnt);
    nCoeff = size(Y,1);
end

%
% THE PREPHEREING STEP
%
if opts.covariance
    
    % calculate the covariance of the noise ...
    if opts.adaptneighbours 
        W = extractNeighbourhood(band1n,opts.bigblock,band2n,prnt);
        W = W(:,index);
    else
        W = extractNeighbourhood(band1n,opts.block,band2n,prnt);
    end
    C_w = innerProd(W)/nCoeff;
    
    S = chol(C_w);
    iS = inv(S);
    
    % ... and of the signal ...
    C_y = innerProd(Y)/nCoeff;
    C_x = C_y - C_w;
    
    % correct possible negative EVs. The overall sum of eigenvalues
    % will remain the same
    [Q,L] = eig(C_x);
    L = diag(diag(L).*(diag(L)>0))*sum(diag(L))/(sum(diag(L).*(diag(L)>0))+(sum(diag(L).*(diag(L)>0))==0));
    C_x = Q*L*Q';
    
    % if all the eigenvalues are zero (or even below) the signal is not
    % detecable at all in the subband therefore it is replaced by zeros
    if sum(L(:)) <= 0
        if verbose
            warning('signal is not detectable in noisy subband\n');
        end
        newBand = zeros(size(band1));
        return;
    end

    % ... transform to whitened noise space ...
    [Q,L] = eig(iS' * C_x * iS);

    % ... the double diagonalization M ...
    M = Q'*iS';
    
    % ... now the noise is white M * C_w * M' = eye...
    nVar = 1;
    nVarInv = 1/nVar;
    
    clear W Q L C_w C_y C_x S dd iS;
else
    M = eye(D,D);
    % pgehler : check which variance to take here. 
    % nVar = noiseVar;
    nVar = var2(band1n);
    nVarInv = 1/nVar;
end

% ... transform the signal
Y = real(M*Y');

% Moment estimation starts ...
m2 = zeros(nCoeff,D);
m4 = zeros(nCoeff,D);

global gl_avg2 gl_avg4

YY = reshape(Y,size(Y,1),sqrt(size(Y,2)),sqrt(size(Y,2)));
for n=1:D
    
    % ... the square band ...
    YY2 = squeeze(YY(n,:,:)).^2;
    if size(gl_avg2) ~= size(YY2)
        gl_avg2 = filter2(ones(opts.Nsize2),ones(size(YY2)));
        if opts.Nsize2 == opts.Nsize4
            gl_avg4 =gl_avg2;
        else
            gl_avg4 = filter2(ones(opts.Nsize4),ones(size(YY2)));
        end
    end
    % ... is locally averaged ...
    mom2 = filter2(ones(opts.Nsize2,opts.Nsize2),YY2);

    % ... and the estimate is corrected for the number of occurences ..
    mom2 = mom2./gl_avg2;
    YY2 = YY2.^2;
    
    % ... same for the fourth moment ...
    mom4 = filter2(ones(opts.Nsize4,opts.Nsize4),YY2);
    mom4 = mom4./gl_avg4;
    m2(:,n) = mom2(:);
    m4(:,n) = mom4(:);
end
clear YY mom2 mom4  % free up some mem

% All coefficients where the signal is not detectable (the signal variance 
% is smaller than the noise variance) are removed from consideration and
% are also excluded from the neighbourhoods they would appear in.
sig = max(m2-nVar,0);
sigzero = sig <= eps;
bigD = repmat(D,size(sig,1),1);
bigD = bigD - row_sum(sigzero);
deadzone = bigD == 0;

if verbose > 2
    for i = 1:max(bigD(:));
        fprintf('%d experts using %d neighbours\n',sum(bigD==i),i);
    end
end


% All patches consiting solely of noise are removed. This is the so-called
% deadzone. 
Y(:,deadzone) = [];
m2(deadzone,:) = [];
m4(deadzone,:) = [];
sig(deadzone,:) = [];
bigD(deadzone) = [];


% Now the parameter adaption starts, first alpha ...
if opts.adaptalpha
    alpha = opts.alphadef * ones(1,size(Y,2));
    % ... binary search over the interval ...
    aa = 0.05:0.01:1;
    intervalLength = size(aa,2);

    % ... get the index of the default value ...
    [dummy,tind] = max(aa==opts.alphadef);
    IN = tind * ones(size(Y,2),1); % ... prepare the index variable ...
    bb = 1:size(IN,1);
    term = sum((m4 - 6*m2*nVar +3*nVar^2)./ ((m2-nVar).^2).*(sig~=0),2); % ... precompute this guy ...
    clear m2 m4     % ...  and free up some mem 
    for tD = 1:max(bigD(:))
        ii = bigD == tD; % ... how many experts with 'tD' neighbours ?
        if sum(ii)==0 % ... don't need to do something ...
            continue;
        end
        % ... for memory reasons we use only one variable ...
        df = repmat(tD^2 * (gamma((tD+4)/2./aa) .* gamma(tD/2./aa))./ (gamma((tD+2)/2./aa).^2),sum(ii),1);

        % ... we approximate with D*(D-1) here instead of compting the nasty
        % correct term ...
        df = df-repmat(term(ii) + tD*(tD-1),1,intervalLength);
        ind = (df(:,1) >0) & (df(:,end)<0);	% ... check if zero is reached ...
        cc = bb(ii);
        [m,IN(cc(ind))] = min(abs(df(ind,:)),[],2); % ... and find the values at this points ...
    end
    % ... now set the alpha to the value found above ...
    alpha = aa(IN);  
    clear m IN term1 term2 term ind;
else
    % no alpha adaption
    alpha = opts.alphadef;
end

% ... estimation of W ... 
warning off MATLAB:divideByZero
term1 = gamma((bigD+2)./(2*alpha'))./gamma(bigD./(2*alpha'))./bigD;
term1 = repmat(term1,1,size(sig,2));
W = term1./sig;
clear term1;

% ... whereever W is infite we set it to zero and the corresponding wavelet
% coefficients as well. One has to include them in the algorithm because
% the denoising takes place in the sphered domain. When unsphereing the
% other coefficients can make a contribution to the center even though in
% the sphered domain its signal is not detectable. 
Winf = isinf(W);
W(Winf) = 0;
z = Y';
clear Y;

% ... precompute some auxiliary variables ...
W2 = 2*W;
DD = ones(D,1);
Kv = nVarInv * z;
alpha1 = alpha'-1;

for t = 1:opts.mapsteps
    z(Winf) = 0;
    gam = alpha' .* (W.*(z.^2) * DD).^(alpha1);
    z = Kv./(W2.*repmat(gam,1,D) + nVarInv);
end
z(Winf) = 0;
warning on MATLAB:divideByZero

% ... sphere back to the original domain ...
z = M\(z');

% ... and    retain the center coefficient only 
Ynew = zeros(nCoeff,1);
Ynew(~deadzone) = z(center,:); 

% ... the borders of the image are not being denoised, it was extended
% beforehand so this has no impact ...
newBand = band1;
newBand(1+Ly:nblv+Ly,1+Lx:nblh+Lx) = reshape(real(Ynew),nblv,nblh);


function s = row_sum(x)
% ROW_SUM   Sum for each row.
% A faster and more readable alternative to sum(x,2).
s = x*ones(size(x,2),1);