explore(fun,x,y)

%Interactive function explorer.
%
%Synopsis:
%
%  explore(fun)
%     Explore function z=fun(x,y) by displaying the result z
%     as a scaled grayscale image. Use the default |x,y| < 2.
%
%  explore(fun,x,y)
%     Start with user-defined coordinates (x,y). Either two
%     vectors or two matrices in the NDGRID format.
%
%  explore
%     Bring last explorer figure on top.
%
%User interface commands (keys):
% 
%  ESC:  Cancel action.
%  E,F:  Export to file.
%  M,P:  Move/pan region.
%  R,T:  Rotate/turn region.
%  S:    Sub-sample image.
%  W:    Display coordinate.
%  Z:    Zoom in/out.
%  ?:    Print help.

%        Marcel Leutenegger  January 2008
%
function explore(fun,x,y)
o=get(0,'CallbackObject');
if isempty(o)
   if nargin
      if ischar(fun)
         fun=eval(['@' fun]);
      end
      if nargin < 2
         [x,y]=ndgrid(-2:0.01:2,-2:0.01:2);
      elseif ~isequal(size(x),size(y)) | any([size(x) size(y)] == 1)
         [x,y]=ndgrid(x,y);
      end
      f=create(fun,x,y);
   else
      s='ShowHiddenHandles';
      t=get(0,s);
      set(0,s,'on');
      f=findobj(get(0,'Children'),'flat','Tag','FunctionExplorer')
      set(0,s,t);
      figure(f(1));
   end
else
   f=o;
   while get(f,'Parent') ~= 0
      f=get(f,'Parent');
   end
   try
      if o == f
         feval(fun,f);
      else
         feval(get(o,'Tag'),f,o);
      end
   catch
      fprintf('%s\n',lasterr);
   end
end


%Create the user interface.
%
% fun    Function handle
% x      Coordinates for
% y      function plot
%
function f=create(fun,x,y)
[m,n]=size(x);
s=get(0,'ScreenSize');
s=s(3:4) - [m n+17];
u.do='No';
u.now='move';
u.next=u.now;
u.fun=fun;
u.dx=step(x);
u.dy=step(y);
u.x=x;
u.y=y;
%
% Figure, axes and image
%
f=figure('Colormap',gray(256),'DoubleBuffer','on','HandleVisibility','off','KeyPressFcn','explore(''key'')','Menubar','none','Name','Function explorer','NumberTitle','off','Position',[s/2 m n+18],'ResizeFcn','explore(''resize'')','Tag','FunctionExplorer','WindowButtonUpFcn','explore(''done'')');
a=axes('CLim',[0 255],'DataAspectRatio',[1 1 1],'Parent',f,'Position',[0 0 1 1]);
i=image(zeros(m,n),'ButtonDownFcn','explore','CDataMapping','scaled','Parent',a,'Tag','init');
set(a,'ButtonDownFcn','explore','Tag','init','YDir','normal','Visible','off');
draw(f,a,u);


%Terminate mouse drag.
%
% f      Figure handle
%
function done(f)
set(f,'Pointer','arrow','WindowButtonMotionFcn','');
a=get(f,'CurrentAxes');
u=get(f,'UserData');
switch u.now
case 'cancel'
   u.now=u.next;
   [x,y]=size(u.x);
   set(a,'XLim',0.5+[0 x],'YLim',0.5+[0 y]);
   set(f,'UserData',u);
   return;
case 'move'
   [x,y]=move(f,u);
   u.x=u.x + [x y]*u.dx;
   u.y=u.y + [x y]*u.dy;
case 'rotate'
   t=turn(f,u);
   [xc,yc]=position(u,u.xc,u.yc);
   x=u.x - xc;
   y=u.y - yc;
   u.x=xc + cos(t)*x + sin(t)*y;
   u.y=yc - sin(t)*x + cos(t)*y;
   u.dx=step(u.x);
   u.dy=step(u.y);
case 'zoom'
   z=zoom(f,u);
   [x,y]=position(u,u.xc,u.yc);
   u.x=x + z*(u.x - x);
   u.y=y + z*(u.y - y);
   u.dx=z*u.dx;
   u.dy=z*u.dy;
end
draw(f,a,u);


%Drag the mouse cursor.
%
% f      Figure handle
%
function drag(f)
a=get(f,'CurrentAxes');
u=get(f,'UserData');
switch u.now
case 'move'
   [x,y]=move(f,u);
   [m,n]=size(u.x);
   x=0.5 + x + [0 m];
   y=0.5 + y + [0 n];
   set(a,'XLim',x,'YLim',y);
case 'rotate'
   t=turn(f,u);
   h=get(a,'Children');
   c=get(h,'UserData');
   [m,n]=size(u.x);
   [X,Y]=ndgrid((1:m)-u.xc,(1:n)-u.yc);
   x=X*cos(t) + Y*sin(t);
   y=Y*cos(t) - X*sin(t);
   for m=1:size(c,3)
      c(:,:,m)=interpn(X,Y,c(:,:,m),x,y);
   end
   if m > 1
      c=uint8(255/max(c(:))*c);
   end
   set(h,'CData',permute(c,[2 1 3]));
case 'zoom'
   z=zoom(f,u);
   [m,n]=size(u.x);
   x=u.xc + z*(0.5 - u.xc + [0 m]);
   y=u.yc + z*(0.5 - u.yc + [0 n]);
   set(a,'XLim',x,'YLim',y);
end


%Calculate and (re)draw image.
%
% f      Figure handle
% a      Axes handle
% u      User data
%
function draw(f,a,u)
set(f,'Pointer','watch');
c=feval(u.fun,u.x,u.y);
h=get(a,'Children');
[x,y]=size(u.x);
if ~isreal(c)
   c=abs(c);
end
c=permute(c,[2 1 3]);
set(a,'XLim',0.5+[0 x],'YLim',0.5+[0 y]);
set(h,'CData',c,'XData',[1 x],'YData',[1 y],'UserData',[]);
set(f,'Pointer','arrow','UserData',u);
drawnow;


%Initialize mouse drag.
%
% f      Figure handle
% o      Image handle
%
function init(f,o)
u=get(f,'UserData');
a=get(f,'CurrentAxes');
p=get(a,'CurrentPoint');
u.now=u.next;
u.point=get(f,'CurrentPoint');
u.xc=p(1);
u.yc=p(3);
switch u.now
case 'move'
   p='cross';
case 'rotate'
   p='left';
   h=get(a,'Children');
   set(h,'UserData',double(permute(get(h,'CData'),[2 1 3])));
case 'zoom'
   p='top';
otherwise
   p='crosshair';
end
set(f,'Pointer',p,'WindowButtonMotionFcn','explore(''drag'')','UserData',u);


%Select the next mouse action.
%
% f      Figure handle
%
function key(f)
u=get(f,'UserData');
switch get(f,'CurrentCharacter')
case 27                    % undo mouse drag
   if ~isempty(get(f,'WindowButtonMotionFcn'))
      u.now='cancel';
      set(f,'UserData',u);
      done(f);
   end
case '?'
   help('explore');
case {'e','E','f','F'}     % export
   here=cd;
   if isfield(u,'path')
      cd(u.path);
   end
   [file,path]=uiputfile({'*.png','Portable Network Graphics (*.png)'; ...
      '*.tif;*.tiff','Tagged Image File Format (*.tif;*.tiff)'; ...
      '*.jpg;*.jpeg','JPEG compliant (*.jpg;*.jpeg)';
      '*.*','All files'},'Save image as');
   cd(here);
   if ~isequal(file,0)
      u.path=path;
      h=get(get(f,'CurrentAxes'),'Children');
      imwrite(get(h,'CData'),fullfile(path,file));
   end
case {'m','M','p','P'}     % move/pan
   u.next='move';
   fprintf('Move/pan.\n');
case {'r','R','t','T'}     % rotate/turn
   u.next='rotate';
   fprintf('Rotate/turn.\n');
case {'s','S'}             % sub-sample
   if isempty(get(f,'WindowButtonMotionFcn'))
      if ~isequal(u.do,'Always')
         u.do=questdlg({'This will apply a 16x16 sub-sampling to the image.','Depending upon the function speed, it may take','a long time to finish.',' ','Do you really want do proceed?'},'Sub-sample image','Yes','No','Always',u.do);
      end
      if ~isequal(u.do,'No')
         subsample(f,u);
      end
   end
case {'w','W'}             % coordinate
   t=get(0,'Unit');
   set(0,'Unit','pixel');
   p=get(0,'PointerLocation');
   set(0,'Unit',t);
   s=get(f,'Position');
   [x,y]=position(u,p(1)-s(1),p(2)-s(2));
   fprintf('(x,y) = (%f,%f)\n',x,y);
case {'z','Z'}             % zoom in/out
   u.next='zoom';
   fprintf('Zoom in/out.\n');
end
set(f,'UserData',u);


%Pan distance.
%
% f      Figure handle
% u      User data
% x,y    Pan distance
%
function [x,y]=move(f,u)
y=u.point - get(f,'CurrentPoint');
x=y(1);
y=y(2);


%Translate mouse position into function coordinate.
%
% xc,yc  Axis coordinate (mouse pointer)
% u      User data
% x,y    Function coordinate
%
function [x,y]=position(u,xc,yc)
x=interpn(u.x,xc,yc);
y=interpn(u.y,xc,yc);


%Resize window - enlarge region.
%
% f      Figure handle
%
function resize(f)
a=get(f,'CurrentAxes');
p=get(f,'Position');
u=get(f,'UserData');
m=max(1,p(3)-1);
n=max(1,p(4)-1);
m=0.5-m/2:m/2;
n=0.5-n/2:n/2;
[x,y]=ndgrid(mean(u.x(:)) + u.dx(1)*m,u.dx(2)*n);
u.x=x + y;
[x,y]=ndgrid(mean(u.y(:)) + u.dy(1)*m,u.dy(2)*n);
u.y=x + y;
draw(f,a,u);


%Sub-sampling positions.
%
% D      Sample step
% d      Sub-sampling
%
function d=sample(D)
[x,y]=ndgrid((-7.5:8)/16,(-7.5:8)/16);
d=[x(:) y(:)]*D;


%Coordinate increment.
%
% r      Coordinate matrix
% d      Increment in both dimensions
%
function d=step(r)
d=[r(2);r(1,2)] - r(1);


%Sub-sample the current region.
%
% f      Figure handle
% u      User data
%
function subsample(f,u)
h=get(get(f,'CurrentAxes'),'Children');
c=get(h,'CData');
try
   set(f,'Interruptible','off','Pointer','watch');
   fprintf('Sub-sampling image (Ctrl-C to escape) ... ');
   dx=sample(u.dx);
   dy=sample(u.dy);
   [m,n]=size(u.x);
   t=time;
   for n=1:n
      [x,y]=ndgrid(dx,u.x(:,n));
      x=x + y;
      [y,z]=ndgrid(dy,u.y(:,n));
      y=y + z;
      z=feval(u.fun,x,y);
      c(n,:,:)=sum(z,1)/m;
      T=time;
      if T - t > 1
         set(h,'CData',c);
         drawnow;
         t=T;
      end
   end
   fprintf('done.\n');
catch
   set(h,'CData',c);
   fprintf('quit.\n');
end
set(f,'Interruptible','on','Pointer','arrow');
drawnow;


%Current time.
%
% t      Time [s]
%
function t=time
t=86400*now;


%Rotation angle.
%
% f      Figure handle
% u      User data
% t      Rotation angle
%
function t=turn(f,u)
t=u.point - get(f,'CurrentPoint');
t=pi/180*t(1);


%Zoom factor.
%
% f      Figure handle
% u      User data
% z      Zoom factor
%
function z=zoom(f,u)
z=get(f,'CurrentPoint') - u.point;
z=pow2(z(2)/50);