function demo_bspline_bivariate
% graphical input of knots 
% for a bivariate B-spline 

global AMB_Data;
AMB_Data=[];
AMB_Data.tau1=[];
AMB_Data.tau2=[];
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.axis=[0,0,0];
AMB_Data.uicontrols=[];
AMB_Data.rotate=0;
AMB_Data.handles=[0];
clf reset;
set(gcf,'WindowButtonDownFcn',@myDownfct);
set(gcf,'WindowButtonMotionFcn',@myMovefct);
set(gcf,'WindowButtonUpFcn',@myUpfct);

% Text for title
uicontrol('Units','normalized', ...
	  'BackgroundColor',get(gcf,'Color'), ...
	  'HorizontalAlignment','center', ...
	  'FontWeight','bold',...
	  'FontSize',20,...
	  'Position',[.05,.92,.9,.07],...
	  'String','Bivariate B-Spline',...
	  'Style','text');

% Text for mouse
AMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'ForegroundColor',[0,0,0], ...
        'FontWeight','bold',...
        'HorizontalAlignment','right', ...
	'Position',[.05,0,.9,.05],...
	'String','Mouse: left = add, right = move, middle = delete knot ',...
	'Style','text');

% Button for new B-spline
uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Restart', ...
	'Position',[.75,.87,.2,.05], ...
	'Callback',@restart,...
	'Style','pushbutton');

% Button for rotate
AMB_Data.uicontrols(2)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Rotate', ...
	'Position',[.05,.87,.2,.05], ...
	'Callback',@rotate,...
	'Style','togglebutton');

% Axis for B-spline
AMB_Data.axis(1)=subplot('position',[.15,.15,.8,.7]);
axis([0,1,0,1,0,1]);
view(3);
axis off;
colormap(winter)
hold on;

% Horizontal axis for knots
AMB_Data.axis(2)=subplot('position',[.15,.06,.8,.05]);
plot([0 1],[0,0],'-k');
axis([-.01,1.01,-.1,.1]);
axis manual
axis off

% Vertical axis for knots
AMB_Data.axis(3)=subplot('position',[.06,.15,.05,.7]);
plot([0 0],[0,1],'-k');
axis([-.1,.1,-.01,1.01]);
axis manual
axis off

%%% Callback Functions %%%

function myDownfct(varargin)
global AMB_Data;

if strcmp(get(gcf,'SelectionType'),'normal') % left mouse 
  p=get(AMB_Data.axis(2),'CurrentPoint');
  p=p(1,1:2);  
  if p(1) >= 0 & p(1) <=1 & p(2)>-.1 & p(2) < .1 
    AMB_Data.tau1=[AMB_Data.tau1;p(1)];
    [AMB_Data.tau1,sind]=sort(AMB_Data.tau1);
    AMB_Data.moveind=find(sind==size(AMB_Data.tau1,1));	  
    AMB_Data.movetype=1;
    myplot;
  else 
    p=get(AMB_Data.axis(3),'CurrentPoint');
    p=p(1,1:2);  
    if p(2) >= 0 & p(2) <=1 & p(1)>-.1 & p(1) < .1 
      AMB_Data.tau2=[AMB_Data.tau2;p(2)];
      [AMB_Data.tau2,sind]=sort(AMB_Data.tau2);
      AMB_Data.moveind=find(sind==size(AMB_Data.tau2,1));	  
      AMB_Data.movetype=2;
    end 
    myplot;	
  end
end

if strcmp(get(gcf,'SelectionType'),'alt') % right mouse 
  p=get(AMB_Data.axis(2),'CurrentPoint');
  p=p(1,1:2);  
  if p(1) >= 0 & p(1) <=1 & p(2)>-.1 & p(2) < .1 
    if ~isempty(AMB_Data.tau1) 
      % nearest knot
      [m,ind]=min(abs(p(1)-AMB_Data.tau1));
      if m < .05 
	AMB_Data.movetype=1; 
	AMB_Data.tau1(ind)=p(1);
	[AMB_Data.tau1,sind]=sort(AMB_Data.tau1);
	AMB_Data.moveind=sind(ind);	  
	myplot;      
      end      	
    end 
  else 
    p=get(AMB_Data.axis(3),'CurrentPoint');
    p=p(1,1:2);  
    if p(2) >= 0 & p(2) <=1 & p(1)>-.1 & p(1) < .1 
      if ~isempty(AMB_Data.tau2) 
	[m,ind]=min(abs(p(2)-AMB_Data.tau2));
	if m < .05 
	  AMB_Data.movetype=2; 
	  AMB_Data.tau2(ind)=p(2);
	  [AMB_Data.tau2,sind]=sort(AMB_Data.tau2);
	  AMB_Data.moveind=sind(ind);	  
	  myplot;      
	end      	
      end 
    end
  end
end

if strcmp(get(gcf,'SelectionType'),'extend') % middle mouse 
  p=get(AMB_Data.axis(2),'CurrentPoint');
  p=p(1,1:2);  
  if p(1) >= 0 & p(1) <=1 & p(2)>-.1 & p(2) < .1 
    if ~isempty(AMB_Data.tau1) 
      % nearest knot 
      [m,ind]=min(abs(p(1)-AMB_Data.tau1));
      if m < .05 
	AMB_Data.tau1(ind)=[];
	myplot;      
      end      	
    end  
  else
    p=get(AMB_Data.axis(3),'CurrentPoint');
    p=p(1,1:2);  
    if p(2) >= 0 & p(2) <=1 & p(1)>-.1 & p(1) < .1 
      if ~isempty(AMB_Data.tau2) 
	% nearest knot 
	[m,ind]=min(abs(p(2)-AMB_Data.tau2));
	if m < .05 
	  AMB_Data.tau2(ind)=[];
	  myplot;   
	end      	
      end  
    end
  end
end
  
function myUpfct(varargin)
global AMB_Data;
AMB_Data.movetype=0;
AMB_Data.moveind=0;
myplot;      

function myMovefct(varargin)
global AMB_Data;
if AMB_Data.movetype==1
  p=get(AMB_Data.axis(2),'CurrentPoint');
  p=p(1,1);
  ind = AMB_Data.moveind;
  tau=[0;AMB_Data.tau1;1];
  if p < tau(ind), p=tau(ind); end
  if p > tau(ind+2), p=tau(ind+2); end 
  AMB_Data.tau1(ind)=p;
  myplot;  
elseif  AMB_Data.movetype==2
  p=get(AMB_Data.axis(3),'CurrentPoint');
  p=p(1,2);
  ind = AMB_Data.moveind;
  tau=[0;AMB_Data.tau2;1];
  if p < tau(ind), p=tau(ind); end
  if p > tau(ind+2), p=tau(ind+2); end 
  AMB_Data.tau2(ind)=p;
  myplot;
end     

function rotate(varargin)
global AMB_Data;
if AMB_Data.rotate
  AMB_Data.rotate=0;
  rotate3d off;
  set(AMB_Data.uicontrols(1),'String','Mouse: left = add, right = move, middle = delete knot ');
else
  AMB_Data.rotate=1;
  rotate3d on;
  set(AMB_Data.uicontrols(1),'String','disable rotation to change knots ');
end

function restart(varargin)
global AMB_Data;
AMB_Data.tau1=[];
AMB_Data.tau2=[];
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.rotate=0;
rotate3d off;
set(AMB_Data.uicontrols(1),'String','Mouse: left = add, right = move, middle = delete knot ');
set(AMB_Data.uicontrols(2),'Value',0);
myplot;

%%% Internal Functions %%%

function myplot
global AMB_Data;
tau1=AMB_Data.tau1;

% calculate multiplicities
dtau1=diff(tau1);
ind1=find(dtau1<1000*eps);
mult=ones(size(tau1));
for k=fliplr(ind1')
  mult(k)=mult(k)+mult(k+1);
  mult(k+1)=[];
  tau1(k+1)=[];
end

% plot horizontal knots
subplot(AMB_Data.axis(2));
cla;
view(2);
plot([0 1],[0,0],'-k');
axis manual
axis([-.01,1.01,-.1,.1]);
axis off
hold on;
n1=size(AMB_Data.tau1,1);
for k=1:length(tau1) 
  plot(tau1(k)*[1 1],[-.08,.08],'-k');
  if mult(k)>1
    text(tau1(k)-.002,.12,sprintf('%d',mult(k)));
  end
end    
if AMB_Data.movetype ==1 
  plot(AMB_Data.tau1(AMB_Data.moveind)*[1 1],[-.08,.08],'-r');
end  

tau2=AMB_Data.tau2;

% calculate multiplicities
dtau2=diff(tau2);
ind2=find(dtau2<1000*eps);
mult=ones(size(tau2));
for k=fliplr(ind2')
  mult(k)=mult(k)+mult(k+1);
  mult(k+1)=[];
  tau2(k+1)=[];
end

% plot vertical knots 
subplot(AMB_Data.axis(3));
cla;
view(2);
plot([0 0],[0,1],'-k');
axis manual
axis([-.1,.1,-.01,1.01]);
axis off
hold on;
n2=size(AMB_Data.tau2,1);
for k=1:length(tau2) 
  plot([-.08,.08],tau2(k)*[1 1],'-k');
  if mult(k)>1
    text(.12,tau2(k)-.002,sprintf('%d',mult(k)));
  end
end    
if AMB_Data.movetype ==2 
  plot([-.08,.08],AMB_Data.tau2(AMB_Data.moveind)*[1 1],'-r');
end  

% plot B-spline
subplot(AMB_Data.axis(1));
cla;
axis([0,1,0,1,0,1])
axis off 
tau1=AMB_Data.tau1;
tau2=AMB_Data.tau2;
if n1 > 1 & n2 > 1 
  t=linspace(0,1,301);
  b1=bspline(tau1,t);
  ind1=find(abs(diff(b1)) > .95);
  b1(ind1)=NaN;
  b2=bspline(tau2,t);
  ind2=find(abs(diff(b2)) > .95);
  b2(ind2)=NaN;
  surf(t,t,b2'*b1,'edgecolor','none','facecolor','interp','facealpha',.5);
  hold on; 
  mesh(tau1,tau2,0*tau2*tau1','facecolor','none','edgecolor','k');
end

%%% Functions from AMB_Programs (included for standalone use) %%%

function b = bspline(xi,x) 
% BSPLINE: evaluation of a B-spline  
% b = bspline(xi,x) 
% values b(j) of a B-spline with knots xi at x(j)

b = 0*x; 
if length(xi) == 2
   % constant B-spline
   b( x>=xi(1) & x<xi(2) ) = 1;
else
   % recursion, discarding terms with zero denominator
   d = xi(end-1)-xi(1);
   if d > 0
      b = b + (x-xi(1)).*bspline(xi(1:end-1),x)/d;
   end
   d = xi(end)-xi(2);
   if d > 0
      b = b + (xi(end)-x).*bspline(xi(2:end),x)/d;
   end
end

