% AGC Automatic Gain Control for geophysical images
%
% Usage: agcim = agc(im, sigma, p, r)
%
% Arguments: im - The input image. NaNs in the image are handled
% automatically.
% sigma - The standard deviation of the Gaussian filter used to
% determine local image mean values and to perform the
% summation used to determine the local gain values.
% Sigma is specified in pixels, try experimenting with a
% wide range of values.
% p - The power used to compute the p-norm of the local image
% region. The gain is obtained from the p-norm. If
% unspecified its value defaults to 2 and r defaults to 0.5
% r - Normally r = 1/p (and it defaults to this) but it can
% specified separately to achieve specific results. See
% Rajagopalan's papers.
%
% Returns: agcim - The Automatic Gain Controlled image.
%
% The algorithm is based on Shanti Rajagopalan's papers, referenced below, with
% a couple of differences.
%
% 1) The gain is computed from the difference between the image and its local
% mean. The aim of this is to avoid any local base-level issues and to allow
% the code to be applied to a wide range of image types, not just magnetic
% gradient data. The gain is applied to the difference between the image and
% its local mean to obtain the final AGC image.
%
% 2) The computation of the local mean and the summation operations used to
% compute the local gain is performed using Gaussian smoothing. The aim of this
% is to avoid abrupt changes in gain as the summation window is moved across the
% image. The effective window size is controlled by the value of sigma used to
% specify the Gaussian filter.
%
% References:
% * Shanti Rajagopalan "The use of 'Automatic Gain Control' to Display Vertical
% Magnetic Gradient Data". 5th ASEG Conference 1987. pp 166-169
%
% * Shanti Rajagopalan and Peter Milligan. "Image Enhancement of Aeromagnetic Data
% using Automatic Gain Control". Exploration Geophysics (1995) 25. pp 173-178
% Copyright (c) 2012 Peter Kovesi
% Centre for Exploration Targeting
% The University of Western Australia
% peter.kovesi at uwa edu au
%
% 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.
% April 2012 - Original version
function agcim = agc(im, sigma, p, r)
% Default values for p and r
if ~exist('p', 'var'), p = 2; end
if exist('p','var') & ~exist('r', 'var')
r = 1/p;
elseif ~exist('r', 'var')
r = 0.5;
end
% Make provision for the image containing NaNs
mask = ~isnan(im);
if any(mask)
im = fillnan(im);
end
% Get local mean by smoothing the image with a Gaussian filter
h = fspecial('gaussian', 6*sigma, sigma);
localMean = filter2(h, im);
% Subtract image from local mean, raise to power 'p' then apply Gaussian
% smoothing filter to obtain a local weighted sum. Finally raise the result
% to power 'r' to obtain the 'gain'. Typically p = 2 and r = 0.5 which will
% make gain equal to the local RMS. The abs() function is used to allow
% for arbitrary 'p' and 'r'.
gain = (filter2(h, abs(im-localMean).^p)).^r;
% Apply inverse gain to the difference between the image and the local
% mean to obtain the final AGC image.
agcim = (im-localMean)./(gain + eps) .* mask;