function [ImageHdl, INTERPOLATED, xpix, ypix] = imagesm( IMAGE, XLIM, YLIM, CRANGE, TYPE, AxesHdl, varargin ) 
% [IMAGEHDL, INTERPOLATED, xpix, ypix] = imagesm( IMAGE, XLIM, YLIM, CRANGE, TYPE, AxesHdl, IMAGE OPTIONS... )
%
% Display an image smoothly interpolated to the resolution of your video display.
% The image is zoomed in to the axis limits specified by XLIM [1x2] and YLIM [1x2],
% and then the selected portion of the image is interpolated 
% using the specified TYPE of interpolation.
%
% Variables:
%
% IMAGEHDL -      Handle of the plotted image object
% INTERPOLATED -  Interpolated grayscale image matrix, before clipping.
%
% IMAGE [NxM] -   Original grayscale image matrix, type double.
% XLIM (default: 0.5+[0 M]) - X Axis limits of the new, interpolated image
% YLIM (default: 0.5+[0 N]) - Y Axis limits of the new, interpolated image
% CRANGE (default: [0 1]) -   Range of useful intensities in IMAGE
% TYPE (default: 'sinc') -    type of interpolation ('sinc', 'linear', 'gaussian', or 'spline')
% AxesHdl (default: gca) -    handle of the axis object on which to plot
% IMAGE OPTIONS - additional options to pass to the matlab image function
%
% An empty matrix ([]) may be given for any of the above (except IMAGE), and defaults will be used.
%
% Zooming algorithm:  The portion of IMAGE selected by XLIM and YLIM is linearly
% interpolated to create a matrix with a number of rows equal to the height in pixels
% of the current axis, and a number of columns equal to the axis width in pixels.
%
% Color normalization algorithm:  IMAGE is normalized so that pixels with a value
% of CRANGE(1) are mapped to the first entry in the current colormap, and pixels 
% with a value of CRANGE(2) or more are mapped to the last entry in the current
% colormap.
%
% Example: Plot the entire image in the matrix IMAGE, mapping values of 255 or
% more to the last entry in the current colormap:
%
% >> imagesm( IMAGE, [], [], [0 255] )
%
% Version 1.5, Mark Hasegawa-Johnson 3/1999.
%

[M, N] = size(IMAGE);

% Figure out what values of XLIM, YLIM, and CRANGE to use
if nargin < 6, AxesHdl = []; end;
if nargin < 5, TYPE = ''; end;
if nargin < 4, CRANGE = []; end
if nargin < 3, YLIM = []; end
if nargin < 2, XLIM = []; end
if length(AxesHdl)==0, AxesHdl=gca; end
if length(TYPE)==0, TYPE='sinc'; end
if length(CRANGE)==0, CRANGE = [0 1]; end
if length(YLIM)==0, YLIM = [ 0.5 M+0.5 ]; end
if length(XLIM)==0, XLIM = [ 0.5 N+0.5 ]; end

% Find the size of the current axes in pixels, and round it to an integer number of pixels.
% Then reset 'Units' to 'normalized', so that the axes will resize when the image does
set(AxesHdl,'Units','pixels');
Hposition = get(AxesHdl,'position');
Hposition(3:4) = round(Hposition(3:4));
set(AxesHdl,'position',Hposition);
set(AxesHdl,'Units','normalized');

% Clip out the desired part of the image
% Add a pixel on every edge, to help valid interpolation
mins = max(1, floor([YLIM(1) XLIM(1)]) - 1);
maxs = ceil([ min(M, YLIM(2)+10) min(N, XLIM(2)+1) ]);
img = IMAGE(mins(1):maxs(1),mins(2):maxs(2));

% Call imresamp to resample the image
[m,n] = size(img);
LX = Hposition(3)/(diff(XLIM)+1);
LY = Hposition(4)/(diff(YLIM)+1);
[ INTERPOLATED, xpix, ypix ] = imresamp( img, LY, LX, TYPE );
% Changed 2/2002 --- find the nearest integer upsampling factor in each direction,
% upsample, then LPF using a first-order Butterworth
%lx=round(LX);
%ly=round(LY);
%xpix=1:(1/lx):n;
%ypix=1:(1/ly):m;
%INTERPOLATED=zeros(ly*m,lx*n);
%INTERPOLATED([0:(m-1)]*ly+1,[0:(n-1)]*lx+1) = lx*ly*img;
%for i=1:m,
%   ay = tan(pi/(2*ly));
%   INTERPOLATED((i-1)*ly+1,:) = filter([1 1],[1+ay, -ay],INTERPOLATED((i-1)*ly+1,:));
%end
%for i=1:(n*lx),
%   ax = tan(pi/(2*lx));
%   INTERPOLATED(:,i) = filter([1 1],[1+ax -ax],INTERPOLATED(:,i));
%end

% Adjust xpix and ypix to match the current image
xpix = xpix - xpix(1) + mins(2);
ypix = ypix - ypix(1) + mins(1);

% Figure out which portion of the interpolated image should be shown
whichx = (xpix >= XLIM(1)) & (xpix <= XLIM(2));
whichy = (ypix >= YLIM(1)) & (ypix <= YLIM(2));
INTERPOLATED = INTERPOLATED(whichy, whichx);
xpix = xpix(whichx);
ypix = ypix(whichy);

% Check to see if the current axes already have an image child
Kids = get(AxesHdl,'Children');
if Kids,
  ImageInd = strmatch('image',get(Kids,'type'));
else,
  ImageInd=[];
end;

% Find out how many colors are in the Colormap, and scale the color range appropriately
NColors = max(size(get(get(AxesHdl,'Parent'),'Colormap')));
INTERPOLATED = NColors * (INTERPOLATED - CRANGE(1))/diff(CRANGE);

% Store previous XLim and YLim in the secret spot searched by matlab zoom('getlimits'), and set axis limits
set(get(AxesHdl,'ZLabel'),'UserData',[get(AxesHdl,'XLim') get(AxesHdl,'YLim')]);
set(AxesHdl, 'XLim', XLIM, 'YLim', YLIM);

% If the current axes already have an image child, add the current CData, XData, and YData
if ImageInd,
  ImageHdl = Kids(ImageInd(1));
  
  % Set new CData
  set(ImageHdl, 'CData', INTERPOLATED, 'XData', xpix, 'YData', ypix, varargin{:});
else,
  % Otherwise, plot a new image on specified axes
  ImageHdl = image(xpix, ypix, INTERPOLATED, 'Parent', AxesHdl, varargin{:});
end;








