Skip to content

@fasiha /texture_filtering_dem.m
Last active

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
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);
@fasiha
Owner

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.