Instantly share code, notes, and snippets.

# fasiha/texture_filtering_dem.m Last active Aug 29, 2015

What would you like to do?
 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)*(xext(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);
Owner

### fasiha commented May 28, 2014

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