www.gusucode.com > 瑞利信道下噪声能量的估计源码程序 > 瑞利信道下噪声能量的估计源码程序/MSP_estimate_version1/ReferenceMethods/Sijbers.m

    % Sijbers's method
% J. Sijbers, et al., Automatic estimation of the noise variance from the histogram of a magnetic resonance image, Physics in medicine and biology, vol. 52.5 (2007): 1335
%
% 06/12/2013
% 
% Tomasz Pieciak
% AGH university of Science and Technology, Krakow, Poland
% pieciak@agh.edu.pl, http://home.agh.edu.pl/pieciak/
%
% ARGUMENTS
%   data - single-coil MRI data
%   bins - histogram bins
%
% FUNCTION RETURNS
%   sigma - estimated noise level (sigma)
%
% USAGE
%    sigma_estimated = Sijbers(data, bins)

function [sigma] = Sijbers(data, bins)

%data = dataset_T1_1mm_noisy;
%window = [5, 5];
%bins = 1000;

% initial noise level searching
[histogram_p, histogram_x] = hist(data(:), bins);
histogram_p_filtered = filtfilt(ones(1,25), 1, histogram_p); % low-pass filter (window 1x25)
[value, index] = max(histogram_p_filtered);
fc = histogram_x(2*index); 

% minimization procedure - least-squares fitting procedure
[n, l] = hist(data(data <= fc), bins); % letters according to eq. 18 in 'Automatic estimation of the noise variance from the histogram of a magnetic resonance image', p. 1340

Nk = sum(n);
K = bins;
F_MIN = @(x) (  Nk*log( exp(-l(1).^2./(2.*x.^2)) - exp(-l(K).^2./(2.*x.^2)) )   -   sum(n(2:K) .* log( exp(-l(1:K-1).^2./(2.*x.^2)) - exp(-l(2:K).^2./(2.*x.^2)) ))  );   
                        % according to eq. 18, p.1340
                        % x - sigma

sigma0 = histogram_x(index); 
[sigma, fval, exitflag] = fminsearch(F_MIN, sigma0);

if(exitflag == 0)
    fprintf('<strong> Sijbers error! </strong>\n')
end