texture_filtering_dem.m

clear();

%% Parameters
% Input file
INFILE = 'ASTGTM2_N43E011_dem.tif';
% Fractional order of the texture filter, between 0 and 2
ALPHA = .5;

%% Helper functions
% Makes a vector of zero-centered Fourier frequency indexes
mkf = @(n) ceil([0:n-1] - n/2);
if numel(mkf) % Test
    assert(all(mkf(3) == [-1 0 1]));
    assert(all(mkf(4) == [-2 -1 0 1]));
end

% Linearly interpolates the min/max of the input to [0, 1]
rescale01 = @(x) interp1(minmax(x(:)'), [0, 1], x);

% Collapses the input array into a column vector
vec = @(x) x(:);

%% The Algorithm
% Read data
t = double(imread(INFILE));

% Compute the size of the FFTs
Nyx = 2.^nextpow2(size(t));
Ny = Nyx(1);
Nx = Nyx(2);

% Make frequency grids and |f|
fxvec = mkf(Nx);
fyvec = mkf(Ny);
[fx fy] = meshgrid(fxvec, fyvec);
absf = sqrt(fx.^2 + fy.^2);

% Spectrum of the image (frequency domain)
tf = fftshift(fft2(t, Ny, Nx)) / numel(t);

% Apply the frequency weighting and convert back to spatial domain
t2 = ifft2(ifftshift(tf .* absf.^ALPHA));

% Get rid of the zero-padding added by power-of-two Nx and Ny
t2 = t2(1:size(t,1), 1:size(t,2));

%% Visualization
% Compute color limits using quantiles of the image AWAY from the edges
edge = 40; % Size of the border to ignore, in pixels
quantiles = [.01/2, .999];
colorlimits = quantile(vec(t2(edge:end-edge-1, edge:end-edge-1))', quantiles);

% Plot the result
figure();
imagesc(t2);
caxis(colorlimits);
axis image;
colormap(gray(256));
colorbar;
title(sprintf('\\alpha=%g', ALPHA));

%% Save the result
% FIXME this is all VERY ugly and requires 
% http://www.mathworks.com/matlabcentral/fileexchange/27959-geotiffwrite

% Helper function that'll clamp values between lower and upper bounds
clip = @(x, ext, val) x.*(x>=ext(1) & x<=ext(2)) + val(1)*(x<ext(1))*val(1) + val(2)*(x>ext(2));

% Clamp the texture to color limits, rescale to 0 to 1000, and convert to int16.
% Here, 0 and 1000 are arbitrary.
outdata = int16(floor(rescale01(clip(t2, colorlimits, colorlimits)) * 1000));

% This is the bounding box from the input GeoTIFF, obtained via gdalinfo.
% FIXME make modular
bbox = [10.9998611,  42.9998611;  12.0001389,  44.0001389] ; %WSEN

% Write a GeoTIF:
% http://www.mathworks.com/matlabcentral/fileexchange/27959-geotiffwrite
opt.bbox = bbox;
geotiffwrite('bar2.tif', bbox, outdata, 16, opt);

Image result and some details about data at http://imgur.com/WgX6qOB.