imshade(x,y,im,angle);

function imshade(x,y,im,angle);
%IMSHADE plot shaded relief with existing colormap
% 
% IMSHADE(X,Y,IM) will plot the matrix IM with x- and y-axis
% vectors X and Y, using the existing colormap with intensity
% scaled to the magnitude of the topographic slope.
%
% IMSHADE(X,Y,IM,ANGLE) will generate shading to approximate
% lighting from the compass direction ANGLE (given in degrees
% clockwise from north).
%
% Plots as a 3-component color image; future changes of colormap
% will not alter image.
%
% Nondirectional shading usually scales well, though it is not very
% robust to spikes in the topographic gradient. The normalization used
% for directional shading does not always yield satisfactory results
% and may benefit from some tweaking.

% Kelsey Jordahl
% Marymount Manhattan College
% Time-stamp: <Fri Oct 30 12:51:13 EDT 2009>

[nrows,ncols]=size(im);
im2=im-min(im(:));
im2=im2*255/max(im2(:));
rgb=ind2rgb(uint8(im2*length(colormap)/255),colormap);
dx=mean(diff(x)); dy=mean(diff(y));
[fx,fy]=gradient(im,dx,dy);
if nargin==4,
  % directional shading
  g=fx*sind(angle)+fy*cosd(angle);
  % might add better normalization step here
  shade=(1-1.5*(g./max(g(:))));                  % flip and scale shading
  shade(shade>1)=1;   shade(shade<0)=0;
else
  % nondirectional shading
  g=sqrt(fx.^2+fy.^2);                    % directionless gradient
  shade=1-(g./max(g(:)));                  % flip and scale shading
end
hsv=rgb2hsv(rgb);
hsv(:,:,3) = shade;                     % set brightness to scaled slope
im3 = hsv2rgb(hsv);
image(x,y,im3);
axis image; axis xy