% NOISECOMP - Function for denoising an image
%
% function cleanimage = noisecomp(image, k, nscale, mult, norient, softness)
%
% Parameters:
% k - No of standard deviations of noise to reject 2-3
% nscale - No of filter scales to use (5-7) - the more scales used
% the more low frequencies are covered
% mult - multiplying factor between scales (2.5-3)
% norient - No of orientations to use (6)
% softness - degree of soft thresholding (0-hard 1-soft)
%
% The convolutions are done via the FFT. Many of the parameters relate
% to the specification of the filters in the frequency plane.
% The parameters are set within the file rather than being specified as
% arguments because they rarely need to be changed - nor are they very
% critical.
%
% Reference:
% Peter Kovesi, "Phase Preserving Denoising of Images".
% The Australian Pattern Recognition Society Conference: DICTA'99.
% December 1999. Perth WA. pp 212-217
% https://www.peterkovesi.com/papers/denoise.pdf
% Copyright (c) 1998-2000 Peter Kovesi
% www.peterkovesi.com
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% September 1998 - original version
% May 1999 -
% May 2000 - modified to allow arbitrary size images
function cleanimage = noisecomp(image, k, nscale, mult, norient, softness)
%nscale = 6; % Number of wavelet scales.
%norient = 6; % Number of filter orientations.
minWaveLength = 2; % Wavelength of smallest scale filter.
%mult = 2; % Scaling factor between successive filters.
sigmaOnf = 0.55; % Ratio of the standard deviation of the Gaussian
% describing the log Gabor filter's transfer function
% in the frequency domain to the filter center frequency.
dThetaOnSigma = 1.; % Ratio of angular interval between filter orientations
% and the standard deviation of the angular Gaussian
% function used to construct filters in the freq. plane.
epsilon = .00001;% Used to prevent division by zero.
thetaSigma = pi/norient/dThetaOnSigma; % Calculate the standard deviation of the
% angular Gaussian function used to
% construct filters in the freq. plane.
imagefft = fft2(image); % Fourier transform of image
[rows,cols] = size(imagefft);
% Create two matrices, x and y. All elements of x have a value equal to its
% x coordinate relative to the centre, elements of y have values equal to
% their y coordinate relative to the centre.
x = ones(rows,1) * (-cols/2 : (cols/2 - 1))/(cols/2);
y = (-rows/2 : (rows/2 - 1))' * ones(1,cols)/(rows/2);
radius = sqrt(x.^2 + y.^2); % Matrix values contain normalised radius from centre.
radius(round(rows/2+1),round(cols/2+1)) = 1; % Get rid of the 0 radius value in the middle so that
% taking the log of the radius will not cause trouble.
theta = atan2(-y,x); % Matrix values contain polar angle.
% (note -ve y is used to give +ve anti-clockwise angles)
clear x; clear y; % save a little memory
sig = [];
estMeanEn = [];
aMean = [];
aSig = [];
totalEnergy = zeros(rows,cols); % response at each orientation.
for o = 1:norient, % For each orientation.
disp(['Processing orientation ' num2str(o)]);
angl = (o-1)*pi/norient; % Calculate filter angle.
wavelength = minWaveLength; % Initialize filter wavelength.
% Pre-compute filter data specific to this orientation
% For each point in the filter matrix calculate the angular distance from the
% specified filter orientation. To overcome the angular wrap-around problem
% sine difference and cosine difference values are first computed and then
% the atan2 function is used to determine angular distance.
ds = sin(theta) * cos(angl) - cos(theta) * sin(angl); % Difference in sine.
dc = cos(theta) * cos(angl) + sin(theta) * sin(angl); % Difference in cosine.
dtheta = abs(atan2(ds,dc)); % Absolute angular distance.
spread = exp((-dtheta.^2) / (2 * thetaSigma^2)); % Calculate the angular filter component.
for s = 1:nscale, % For each scale.
% Construct the filter - first calculate the radial filter component.
fo = 1.0/wavelength; % Centre frequency of filter.
rfo = fo/0.5; % Normalised radius from centre of frequency plane
% corresponding to fo.
logGabor = exp((-(log(radius/rfo)).^2) / (2 * log(sigmaOnf)^2));
logGabor(round(rows/2+1),round(cols/2+1)) = 0; % Set the value at the center of the filter
% back to zero (undo the radius fudge).
filter = logGabor .* spread; % Multiply by the angular spread to get the filter.
filter = fftshift(filter); % Swap quadrants to move zero frequency
% to the corners.
% Convolve image with even an odd filters returning the result in EO
EOfft = imagefft .* filter; % Do the convolution.
EO = ifft2(EOfft); % Back transform.
aEO = abs(EO);
if s == 1
% Estimate the mean and variance in the amplitude response of the smallest scale
% filter pair at this orientation.
% If the noise is Gaussian the amplitude response will have a Rayleigh distribution.
% We calculate the median amplitude response as this is a robust statistic.
% From this we estimate the mean and variance of the Rayleigh distribution
medianEn = median(reshape(aEO,1,rows*cols));
meanEn = medianEn*.5*sqrt(-pi/log(0.5));
RayVar = (4-pi)*(meanEn.^2)/pi;
RayMean = meanEn;
estMeanEn = [estMeanEn meanEn];
sig = [sig sqrt(RayVar)];
%% May want to look at actual distribution on special images
% hist(reshape(aEO,1,rows*cols),100);
% pause(1);
end
% Now apply soft thresholding
T = (RayMean + k*sqrt(RayVar))/(mult^(s-1)); % Noise effect inversely proportional to
% bandwidth/centre frequency.
validEO = aEO > T; % Find where magnitude of energy exceeds noise.
V = softness*T*EO./(aEO + epsilon); % Calculate array of noise vectors to subtract.
V = ~validEO.*EO + validEO.*V; % Adjust noise vectors so that EO values will
% not be negated
EO = EO-V; % Subtract noise vector.
totalEnergy = totalEnergy + EO;
wavelength = wavelength * mult; % Wavelength of next filter
end
end % For each orientation
disp('Estimated mean noise in each orientation')
disp(estMeanEn);
cleanimage = real(totalEnergy);
%imagesc(cleanimage), title('denoised image'), axis image;