Home > Matching_SYM_LSM > src > Functions > Noise_estimation > noise_standard_deviation_estimation.m

noise_standard_deviation_estimation

PURPOSE ^

Estimate noise variance function for each channel of image IMG.

SYNOPSIS ^

function [s,means,sigma] = noise_standard_deviation_estimation(img,K)

DESCRIPTION ^

 Estimate noise variance function for each channel of image IMG.

 img     = image with c channels
 K       = number of variances (default = 10)

 s       = c x 256 [0:1:255] lookup table for noise standard deviation
 means   = Kx1 vector of mean intensities, where variance is stimated
 sigmas  = Kx1 ve ctor of estimated variances

 For each variance to be estimated we select corresponding image pixels
 and robustly estimate their variance from their mean gradient magnitude.
 We iteratively exclude pixels with large gradients that are probably
 caused by image content rather than noise.
 See W. Foerstner, "Image Preprocessing for Feature Extraction in Digital
 Intensity, Color and Range Images" (2000), Eq. (32) for more details.

 For normalizing the image we need to assign a noise variance to each
 discrete gray value. Therefore we perform a spline interpolation of the
 discrete samples from the previous step.

 Taken from noiseNormalize.m
 
 See also: sfop

 Licence:
   For internal use only.

 Authors:
   Timo Dickscheid, Falko Schindler, Wolfgang Foerstner
   Department of Photogrammetry
   Institute of Geodesy and Geoinformation
   University of Bonn
   Bonn, Germany

 Contact person:
   Falko Schindler (falko.schindler@uni-bonn.de)

 Copyright 2012-2012

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [s,means,sigma] = noise_standard_deviation_estimation(img,K)
0002 % Estimate noise variance function for each channel of image IMG.
0003 %
0004 % img     = image with c channels
0005 % K       = number of variances (default = 10)
0006 %
0007 % s       = c x 256 [0:1:255] lookup table for noise standard deviation
0008 % means   = Kx1 vector of mean intensities, where variance is stimated
0009 % sigmas  = Kx1 ve ctor of estimated variances
0010 %
0011 % For each variance to be estimated we select corresponding image pixels
0012 % and robustly estimate their variance from their mean gradient magnitude.
0013 % We iteratively exclude pixels with large gradients that are probably
0014 % caused by image content rather than noise.
0015 % See W. Foerstner, "Image Preprocessing for Feature Extraction in Digital
0016 % Intensity, Color and Range Images" (2000), Eq. (32) for more details.
0017 %
0018 % For normalizing the image we need to assign a noise variance to each
0019 % discrete gray value. Therefore we perform a spline interpolation of the
0020 % discrete samples from the previous step.
0021 %
0022 % Taken from noiseNormalize.m
0023 %
0024 % See also: sfop
0025 %
0026 % Licence:
0027 %   For internal use only.
0028 %
0029 % Authors:
0030 %   Timo Dickscheid, Falko Schindler, Wolfgang Foerstner
0031 %   Department of Photogrammetry
0032 %   Institute of Geodesy and Geoinformation
0033 %   University of Bonn
0034 %   Bonn, Germany
0035 %
0036 % Contact person:
0037 %   Falko Schindler (falko.schindler@uni-bonn.de)
0038 %
0039 % Copyright 2012-2012
0040 
0041 % TODO: How to define number of bins B and discretization (here 256)?
0042 % up to now: input or default = 10
0043 %% default values
0044 K(nargin < 2)       = 10;
0045 
0046 
0047 %% handle multiple channels
0048 if size(img, 3) > 1
0049     means = zeros(size(img, 3), K);
0050     sigma = zeros(size(img, 3), K);
0051     for c = 1 : size(img, 3)
0052         [s(c,:),means(c,:),sigma(c,:)] = noise_standard_deviation_estimation(img(:, :, c),K);
0053     end
0054     return
0055 else
0056     
0057 end
0058 %% check type of image
0059 if max(img(:)) <= 255
0060     % range[0,255]
0061     type_img = 255;
0062     % scale image
0063     img=double(img)/255;
0064 end
0065 if max(img(:)) <= 1
0066     % range [0,1]
0067     type_img = 1;
0068 end
0069 
0070 %% Image gradients
0071 % We will estimate variances from image gradients, thus need to convolve
0072 % the image with two derivative kernels and add the squared gradients.
0073 grad = imfilter(img, [-1, 0, 1], 'replicate').^2 + ...
0074        imfilter(img, [-1; 0; 1], 'replicate').^2;
0075 
0076 %% Variance estimation
0077 % First we sort all pixels and their gradients with ascending gray value.
0078 
0079 [F, idx] = sort(img(:));
0080 
0081 % collect gradients
0082 G = grad(idx);
0083 %%
0084 % We initialize storage for a certain number of variance estimates with
0085 % their corresponding mean gray value.
0086 [means, sigma] = deal(zeros(1, K));
0087 %%
0088 % For each variance to be estimated we select corresponding image pixels
0089 % and robustly estimate their variance from their mean gradient magnitude.
0090 % We iteratively exclude pixels with large gradients that are probably
0091 % caused by image content rather than noise.
0092 % See W. Foerstner, "Image Preprocessing for Feature Extraction in Digital
0093 % Intensity, Color and Range Images" (2000), Eq. (32) for more details.
0094 for b = 1 : numel(means)
0095     selection = round((b - 1) / numel(means) * numel(img) + 1) : ...
0096                 round( b      / numel(means) * numel(img));
0097     g = G(selection);
0098     mu = mean(g);
0099     for iter = 1 : 3
0100         mu = 2.392211191177332 * mean(g(g < mu));
0101     end
0102     sigma(b) = sqrt(mu / 4);
0103     means(b) = mean(F(selection));
0104 end
0105 % regularize means, if image is under- or overexposed
0106 % prior
0107 means0=([1:K]-0.5)/K; w = 0.01;
0108 % weighted mean
0109 means = w*means0 +(1-w)*means;
0110 %% Interpolation
0111 % For normalizing the image we need to assign a noise variance to each
0112 % discrete gray value. Therefore we perform a spline interpolation of the
0113 % discrete samples from the previous step
0114 g_domain = linspace(0, 1, 256);
0115 if K > 1 
0116     s = interp1(means, sigma, g_domain, 'linear');
0117 end;
0118 % and set the values outside the given range to the first and last value
0119 s(g_domain < means(1)) = sigma(1);
0120 s(g_domain > means(end)) = sigma(end);
0121 % limit s from below by rounding error
0122 s=max(sqrt(1/12)/255,s);
0123 % scale s
0124 if type_img == 255
0125     s=s*255;
0126 end
0127 return

Generated on Sun 19-Jul-2020 23:00:25 by m2html © 2005