function demo_bspline
% graphical input of knots and recursion
% for a B-spline 

global AMB_Data;
AMB_Data=[];
AMB_Data.tau=[];
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.axis=[0,0];
AMB_Data.uicontrols=[];
AMB_Data.recursionstep=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','B-Spline',...
	  'Style','text');

% Text for mouse
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 recursion
AMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Show recursion', ...
	'Position',[.05,.87,.2,.05], ...
	'Callback',@recursion,...
	'Style','pushbutton');

% Axis for B-spline
AMB_Data.axis(1)=subplot('position',[.05,.15,.9,.7]);
axis([-.01,1.01,-.03,3.03]);
axis manual;
axis off
hold on;
fill([-.1,1.1,1.1,-.1],[-.1,-.1,3.1,3.1],[1,1,1],'edgecolor','none');
plot([0,1,1,0,0],[0,0,3,3,0],'-k');

% Axis for knots
AMB_Data.axis(2)=subplot('position',[.05,.06,.9,.05]);
plot([0 1],[0,0],'-k');
axis([-.01,1.01,-.1,.1]);
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.tau=[AMB_Data.tau;p(1)];
    [AMB_Data.tau,sind]=sort(AMB_Data.tau);
    AMB_Data.moveind=find(sind==size(AMB_Data.tau,1));	  
    AMB_Data.movetype=1;
    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.tau) 
      % nearest knot
      [m,ind]=min(abs(p(1)-AMB_Data.tau));
      if m < .05 
	AMB_Data.movetype=1; 
	AMB_Data.tau(ind)=p(1);
	[AMB_Data.tau,sind]=sort(AMB_Data.tau);
	AMB_Data.moveind=sind(ind);	  
	myplot;      
      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.tau) 
      % nearest knot 
      [m,ind]=min(abs(p(1)-AMB_Data.tau));
      if m < .05 
	AMB_Data.tau(ind)=[];
	myplot;      
      end      	
    end  
  end
end
  

function myUpfct(varargin)
global AMB_Data;
if AMB_Data.movetype
  AMB_Data.movetype=0;
  AMB_Data.moveind=0;
  myplot;      
end

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

function restart(varargin)
global AMB_Data;
AMB_Data.tau=[];
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.recursionstep=0;
set(AMB_Data.uicontrols(1),'string','Show recursion')
myplot;

function recursion(varargin)
global AMB_Data;
tau=AMB_Data.tau;
n=size(tau,1);

% no recursion possible
if n < 3,return;end


k=2+AMB_Data.recursionstep;
plotspline;

% last step -> end recursion
if (AMB_Data.recursionstep>0) & (n <= k+1)
  AMB_Data.recursionstep=0;    
  set(AMB_Data.uicontrols(1),'string','Show recursion')
  return
end

% plot spline
t=linspace(0,1,501);

% plot lower level
if n > k
  plot([0,1],[1.8,1.8],'-k');
  for l=1:k 
    b=bspline(tau(l:end-k+l),t);
    ind=find(abs(diff(b)) > .95);
    b(ind)=NaN;
    plot(t,1.8+.5*b,'-k');
  end
end
% plot upper level
if n > k+1
  plot([0,1],[2.4,2.4],'-k');
  for l=1:k+1 
    b=bspline(tau(l:end-k-1+l),t);
    ind=find(abs(diff(b)) > .95);
    b(ind)=NaN;
    plot(t,2.4+.5*b,'-k');
  end
end 

AMB_Data.recursionstep=AMB_Data.recursionstep+1;
% update button text
if n > k+2 
  set(AMB_Data.uicontrols(1),'string','Next step')    
else
  set(AMB_Data.uicontrols(1),'string','End recursion')    
end

%%% Internal Functions %%%

function myplot
global AMB_Data;
tau=AMB_Data.tau;

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

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

% plot B-spline
plotspline;
AMB_Data.recursionstep=0;    
set(AMB_Data.uicontrols(1),'string','Show recursion')

function plotspline
global AMB_Data
subplot(AMB_Data.axis(1));
tau=AMB_Data.tau;
cla;
fill([-.1,1.1,1.1,-.1],[-.1,-.1,3.1,3.1],[1,1,1],'edgecolor','none');
plot([0,1,1,0,0],[0,0,3,3,0],'-k');
if length(tau)>1 
  t=linspace(0,1,501);
  b=bspline(tau,t);
  ind=find(abs(diff(b)) > .95);
  b(ind)=NaN;
  plot(t,1.5*b,'-b','linewidth',2);
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


