surfplot(data,wsize,sigma);

function surfplot(data,wsize,sigma);
% SURFPLOT(data,wsize,sigma) - displays a dataset in 3D overlain with a 
% highpass filtered version of itself as colour - best for smooth data
% data = 2D matrix of double
% ksize = kernel size, must be odd, eg 5
% sigma = width of gaussian function derivative in filter eg 3
%
% type surfplot with no arguments to see demo plot
%
% Requires image processing toolbox for fspecial

% ****************************************************************
% *** GRJ Cooper 
% *** School of Geosciences, University of the Witwatersrand
% *** Johannesburg, South Africa
% *** [e-mail]
% *** www.wits.ac.za/science/geophysics/gc.htm
% ****************************************************************
% *** For geophysical datasets, use vertical derivatives as the
% *** high pass filter type
% ****************************************************************

if nargin==0 % load demo data
 load demodata.mat; 
 wsize=3; sigma=1;
end;

[nr,nc]=size(data);
if wsize<3 wsize=3; end; wsize=2*(floor(wsize/2))+1; % Check parameters
if sigma<=0 sigma=wsize/3; end;
w2=floor(wsize/2);

disp('Filtering');
se=fspecial('log',wsize,sigma);
hp=conv2(data,se,'valid');
data1=data(w2+1:nr-w2,w2+1:nc-w2);

% Threshold gradient data

disp('Thresholding');
hpm=mean(hp(:)); hps=std(hp(:));
for I=1:nr-2*w2
 for J=1:nc-2*w2
  if hp(I,J)>hpm+2*hps hp(I,J)=hpm+2*hps; end;
  if hp(I,J)<hpm-2*hps hp(I,J)=hpm-2*hps; end;
 end;
end;

figure(1);
subplot(2,2,1); imagesc(data1); axis equal; axis tight; axis xy; title('Data');
subplot(2,2,2); imagesc(se); axis equal; axis tight; axis xy; title('Filter'); colorbar;
subplot(2,2,3); imagesc(hp); axis equal; axis tight; axis xy; title('Filtered Data'); axis off;
subplot(2,2,4); surf(data); shading flat; axis tight; view([15 85]); title('Data in 3D'); axis off;
colormap('gray');
currfig=get(0,'CurrentFigure'); set(currfig,'numbertitle','off');
set(currfig,'name','Data and Filter');
figure(2);
subplot(1,1,1); surf(data1,hp); shading flat; axis tight; view([15 85]);
title('Data in 3D Overlain with Filtered Data');
colormap('gray'); axis off;
currfig=get(0,'CurrentFigure'); set(currfig,'numbertitle','off');
set(currfig,'name','Enhanced 3D Surface');