Parameters(subb)

function [LB,th] = Parameters(subb)
% Select adaptively the optimal threshold and block size
% using the SURE principle 
% subb: inputted noisy subband
% LB: outputted block size
% th: outputted threshold value
%
% Author: Zhou Dengwen
% [e-mail]
% Department of Computer Science & Technology
% North China Electric Power University(Beijing)(NCEPU)
%
% Reference(s):
% Zhou Dengwen and Shen Xiaoliu, "Image denoising using block
% thresholding," in Proc. 2008 congress on image and signal 
% processing,Sanya, Hainan, ChinaMay 2008, pp. 335C338.
%
% Last time modified: Jun. 11, 2009
% 

% Compute the risk and the corresponding threshold
% of the different block sizes
LBmax = ceil(size(subb,1)^0.75);
p = floor(log2(LBmax));
for i = 1:p
    LB_set(i) = 2^i;
end
k = 0;
for len = LB_set
    k = k + 1;
    [risk(k),thres(k)] = SureBlockThreshold(subb,len);
end

% Obtain the optimal block size and the corresponding threshold
[guess,ibest] = min(risk);
th = thres(ibest); LB = LB_set(ibest);

%---------------------------------------------------------------
function [guess,th] = SureBlockThreshold(subb,LB)
% Select an optimal threshold for a given block size.
% (Note: The threshold lambda in the original paper has
% modified as lambda^2)
% subb: inputted noisy subband
% LB: inputted block size
% guess: outputted minimum risk
% th: outputted threshold value
%

% Compute the square sum of all pixels in a block 
% and generate a threshold series
LB2 = LB*LB;
if LB == 1
    sb = subb(:).^2;
    if min(sb) > 0 
       Thres = [sb; 0];
    end
else
    sb = BlockSquareSum(subb,LB);
    Thres = sb(find(sb >= LB2-2));
    if min(Thres) > LB2-2
       Thres = [Thres; LB2-2];
    end
end

% Compute the threshold corresponding to the minimum risk
risk = zeros(length(Thres),1); m = 0;
blc_num = length(sb); Thres = Thres';

for th = Thres
        
    risk_temp = 0;
    for k = 1:blc_num
        if sb(k) > th
           riskb = LB2+(th*th-2*th*(LB2-2))/sb(k);
        else
           riskb = sb(k)-LB2;
        end
        risk_temp = risk_temp+riskb; 
    end
    m = m + 1; risk(m) = risk_temp;
    
end

[guess,ibest] = min(risk);
th = Thres(ibest);

%---------------------------------------------------------------
function sb = BlockSquareSum(im,L)
% Compute the square sum of the pixels of every block in an image
% L: inputted block size
% im: inputted image
% sb: outputted result
% 

% Compute the number of the blocks
[nRow,nCol] = size(im);
p = mod(nRow, L);
q = mod(nCol, L);
TempIm = padarray(im,[L-p,L-q], 'post'); % filled with 0
blc_row = ceil(nRow/L); 
blc_col = ceil(nCol/L);

% Compute the square sum of the pixels of every block
sb = zeros(blc_row*blc_col,1); m = 0;
for i = 1:blc_row
    for j = 1:blc_col            
        b1 = (i-1)*L; b2 = (j-1)*L;
        bm = TempIm(b1+1:b1+L,b2+1:b2+L);
        m = m + 1; sb(m) = sum(bm(:).^2);               
    end
end