% TILTDERIV Tilt derivative of potential field data
%
% Usage: td = tiltderiv(im, hx, hy, padding)
%
% Arguments: im - Input potential field image.
% hx, hy - Grid element size, defaults to 1, 1
% padding - Width of tapered padding to apply to the image to reduce
% edge effects. Depending on the degree of cyclic
% discontinuity in your data values of up to, say, 100
% can be useful. Defaults to 50.
%
% Returns: td - The tilt derivative.
%
% Reference:
% Hugh G. Miller and Vijay Singh. Potential field tilt - a new concept for
% location of potential field sources. Applied Geophysics (32) 1994. pp
% 213-217.
%
% See also: ANALYTICSIGNAL, FREQDERIV, DEALIAS
% Copyright (c) 2015-2017 Peter Kovesi
% Centre for Exploration Targeting
% The University of Western Australia
% 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.
%
% March 2015
% October 2017 - Changed to use frequency domain horizontal derivatives and
% to allow for data padding. Specify hx, hy.
% December 2017 - Incorporate dealiasing/lowpass filtering after the padding.
% Unwrapped code to avoid repeated FFT operations
function td = tiltderiv(im, hx, hy, padding)
assert(size(im,3) == 1, 'Image must be single channel');
if ~exist('hx', 'var') && ~exist('hy', 'var'), hx = 1; hy = 1; end
if ~exist('padding', 'var'), padding = 50; end
mask = ~isnan(im);
im = impad(fillnan(im), padding, 'cosine');
IM = fft2(im);
% After padding the data it is almost essential to apply a low pass
% filter with a high cut off to avoid artifacts propagating into the
% images. The tilt derivative can be quite 'delicate' in this regard.
%
% The filter is constructed from the product of two low-pass Butterworth
% filters. One filter cuts out high frequency components in the x direction
% and the other in the y direction. The assumption being that image
% artifacts are a result of the gridding used to form the image. Also,
% the padding process, while made as smooth as possible cyclically,
% introduces high frequency content in the padded regions perpendicular
% to the image edges. I believe this causes artifacts in the tilt
% derivative if they are not filtered.
% The two low-pass filters in x and y are defined as follows:
% f1 = 1.0 ./ (1.0 + (ux ./ cx).^(2*n));
% f2 = 1.0 ./ (1.0 + (uy ./ cy).^(2*n));
% f = f1.*f2;
cx = 0.45; cy = 0.45; n = 15; % Cut off at 0.45 and a high order of n
[~, ux, uy] = filtergrid(size(im));
f = 1.0 ./ ((1.0 + (ux ./ cx).^(2*n)) .* (1.0 + (uy ./ cy).^(2*n)));
IM = IM.*f; % Apply the low-pass filter to the image
% Generate wavenumber grids from the filter grids
kx = 2*pi*ux/hx;
ky = 2*pi*uy/hy;
k = sqrt(kx.^2 + ky.^2);
% Take horizontal and vertical first derivatives and remove the padding.
dx = imtrim(real(ifft2(IM .* (i*kx))), padding);
dy = imtrim(real(ifft2(IM .* (i*ky))), padding);
dv = imtrim(real(ifft2(IM .* k)), padding);
% Finally compute the tilt derivative and mask out any NaNs that were in
% the input image.
td = atan(dv./(sqrt(dx.^2 + dy.^2) + eps)) .* double(mask);