% TERRACE - Terracing operator for potential field data.
%
% Useage: terrimg = terrace(img, nIter, sigma)
%
% Arguments:
% img - Input image to be terraced
% nIter - Number of iterations. This controls the degree of
% terracing. Values from 1 to 20 are useful. Small values can
% be used to provide a simple sharpening effect.
% sigma - Standard deviation of Laplacian of Gaussian used to compute the
% 2nd derivatives. This controls the scale of analysis. Large
% values will result in small features being smoothed out. Try
% values from 1 to 5. Typically small values are better.
%
% Returns:
% terrimg - Terraced image.
%
%
% This implementation differs from that suggested by Cordell and McCafferty
% and by Cooper and Cowan in two ways:
%
% 1) The curvature/second derivative information is fixed at the values computed
% from the original image. It is not recomputed on each iteration as this is
% highly unstable and leads to anomalous 'ripples' in the output. Such ripples
% can be seen in figures 4 and 5 of Cooper and Cowan's paper. The key attribute
% of this code is that the output remains stable no matter how many iterations
% are applied.
%
% 2) Curvature is not computed, all that is required is the 2nd derivative. If
% the 2nd derivative is +ve then curvature will be +ve and vice-versa. The
% Laplacian of Gaussian used to compute 2nd derivatives, the standard deviation
% can be used to control the scale of analysis.
% References:
%
% Cordell L. and McCafferty A.E. 1989. A terracing operator for physical
% property mapping with potential field data. Geophysics 54, 621–634.
%
% G.R.J. Cooper and D.R. Cowan ASEG 2009. "Terracing potential field data"
% Copyright (c) 2017 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.
% PK September 2017
function terrimg = terrace(img, nIter, sigma)
% Structuring element for determining local min and max value. (Note a
% large circular structuring element can be used in which case fewer
% iterations are needed. However I think things are more 'stable' with a
% small SE)
se = strel('square',3);
terrimg = img;
% Generate Laplacian of Gaussian filter and apply to image.
sze = ceil(8*sigma);
if ~mod(sze,2), sze = sze+1; end
h = fspecial('log', sze, sigma);
log = filter2(h, terrimg);
% Iterate to sharpen points of maximum gradient and flatten regions into
% terraces.
for i = 1:nIter
% Get min and max value in neighbourhood of each pixel
minv = imerode(terrimg, se);
maxv = imdilate(terrimg, se);
% Replace values at a point of -ve 2nd derivative with the max value in
% the local neighbourhood. Vice-versa for points with +ve 2nd derivative.
ind = find(log < 0);
terrimg(ind) = maxv(ind);
ind = find(log > 0);
terrimg(ind) = minv(ind);
end