function demo_spline_curve
% graphical input of control points and 
% knots for a spline curve 

global AMB_Data;
AMB_Data=[];
AMB_Data.degree=3;
AMB_Data.tau=[0;.1*ones(AMB_Data.degree,1);.9*ones(AMB_Data.degree,1);1];
AMB_Data.c=zeros(0,2);
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.axis=[0,0];
AMB_Data.handles=[0,0,0];
AMB_Data.uicontrols=[];
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','Spline Curve',...
	  'Style','text');

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

% Text for degree
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.05,.85,.1,.05],...
	'String','Degree:',...
	'Style','text');

% Input for degree
AMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','right', ...
	'Position',[.15,.86,.03,.05],...
	'String',sprintf('%1d',AMB_Data.degree),...
	'Style','edit',...
	'Callback',@degree_update);

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

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


%%% Callback Functions %%%
function myDownfct(varargin)
global AMB_Data;

if strcmp(get(gcf,'SelectionType'),'normal') % new
  np=size(AMB_Data.c,1);
  n=AMB_Data.degree;
  p=get(AMB_Data.axis(1),'CurrentPoint'); % curve axis
  p=p(1,1:2);  
  if min(p) >= 0 & max(p) <=1
    AMB_Data.c=[AMB_Data.c;p];    
    % determine new knot
    if  np > n
      f= (np-n)/(np-n+1);
      l=diff(AMB_Data.tau([n+1,end-n])); 
      innertau=(AMB_Data.tau(n+2:end-n-1)-AMB_Data.tau(n+1))/l;       
      innertau=[f*innertau;f];      
      AMB_Data.tau=[AMB_Data.tau(1:n+1);AMB_Data.tau(n+1)+l*innertau;...
	  AMB_Data.tau([end-n:end])];
    end
    AMB_Data.movetype=1;
    AMB_Data.moveind=np+1;
  else
    if np > AMB_Data.degree
      p=get(AMB_Data.axis(2),'CurrentPoint'); % knot axis
      p=p(1,1:2);  
      if p(1) >= AMB_Data.tau(2) & p(1) <= AMB_Data.tau(end-1) & p(2)>-.1 & p(2) < .1 
	[pn,tau]=knot_insertion(AMB_Data.c,AMB_Data.tau,p(1));
	AMB_Data.tau=tau;
	AMB_Data.c=pn;	
	AMB_Data.movetype=2; 
	AMB_Data.moveind=find(AMB_Data.tau==p(1));
      end      	
    end  
  end
  myplot;
end

if strcmp(get(gcf,'SelectionType'),'alt') % move
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);  
  if min(p) >= 0 & max(p) <=1 % move conrol point
    % nearest point
    [m,ind]=min((p(1)-AMB_Data.c(:,1)).^2+(p(2)-AMB_Data.c(:,2)).^2);
    if m < .01 % really near
      AMB_Data.movetype=1; 
      AMB_Data.moveind=ind;
      AMB_Data.c(ind,1:2)=p(1:2);
      myplot;      
    end
  else  
    p=get(AMB_Data.axis(2),'CurrentPoint'); % move knot
    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 % really near
	  AMB_Data.movetype=2; 
	  AMB_Data.tau(ind)=p(1);
	  [AMB_Data.tau,sind]=sort(AMB_Data.tau);
	  AMB_Data.moveind=find(AMB_Data.tau==p(1));
	  myplot;      
	end      	
      end  
    end
  end
end

if strcmp(get(gcf,'SelectionType'),'extend') % middle mouse 
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);  
  if min(p) >= 0 & max(p) <=1 % Middle mouse: delete
    [m,ind]=min((p(1)-AMB_Data.c(:,1)).^2+(p(2)-AMB_Data.c(:,2)).^2);
    if m < .01 
      AMB_Data.c(ind,:)=[];
      n=AMB_Data.degree;
      tau=AMB_Data.tau;      
      if length(tau)>2*n+3
	tau = [tau(1:1+n);(tau(2+n:end-n-2)+tau(3+n:end-n-1))/2;tau(end-n:end)];
      elseif length(tau) == 2*n+3
	tau = [tau(1:1+n);tau(end-n:end)];
      end 
      AMB_Data.tau=tau;
      myplot;      
    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 == 2
  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;
elseif AMB_Data.movetype == 1
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);
  if min(p) >= 0 & max(p) <=1
    ind = AMB_Data.moveind;
    AMB_Data.c(ind,1:2)=p;
    myplot;
  end     
end

function restart(varargin)
global AMB_Data;
AMB_Data.tau=[0;.1*ones(AMB_Data.degree,1);.9*ones(AMB_Data.degree,1);1];
AMB_Data.c=zeros(0,2);
AMB_Data.movetype=0;
AMB_Data.moveind=0;
myplot;

function degree_update(varargin)
global AMB_Data;
degree= str2double(get(AMB_Data.uicontrols(1),'String'));
if degree < 1 | degree > 9
  set(AMB_Data.uicontrols(1),'String',sprintf('%1d',AMB_Data.degree));
else
  AMB_Data.degree=degree;
  restart;
end

%%% Internal Functions %%%
function myplot
global AMB_Data;

n=size(AMB_Data.c,1);

% knot axis
subplot(AMB_Data.axis(2));
cla;
plot([0 1],[0,0],'-k');
axis manual
axis([-.01,1.01,-.1,.1]);
axis off
hold on;
tau=AMB_Data.tau;

% plot D_xi
tau0=tau(AMB_Data.degree+1);
tau1=tau(end-AMB_Data.degree);
fill([tau0,tau0,tau1,tau1],[-.06,.06,.06,-.06],.5*[1,1,1]);

% 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 tau
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  == 2
  plot(AMB_Data.tau(AMB_Data.moveind)*[1 1],[-.08,.08],'-r');
end  

% spline axis
subplot(AMB_Data.axis(1));
np=size(AMB_Data.c,1);
tau=AMB_Data.tau;

% check for regularity
if any(mult > AMB_Data.degree) | ...
      tau(AMB_Data.degree+1) >= tau(AMB_Data.degree+2) | ...
      tau(end - AMB_Data.degree-1) >= tau(end-AMB_Data.degree)
  if AMB_Data.handles(3)==0  
      AMB_Data.handles(3)=...
	  text(.5,.5,sprintf('Knot sequence not %d-regular',AMB_Data.degree),...
	       'horizontalalign','center','fontweight','bold','fontsize',10);
  end
  if AMB_Data.handles(2)~=0  
    delete(AMB_Data.handles(2));
    AMB_Data.handles(2)=0;
  end
  if AMB_Data.handles(1)~=0  
    delete(AMB_Data.handles(1));
    AMB_Data.handles(1)=0;
  end
  return;
else  
  if AMB_Data.handles(3)~=0  
    delete(AMB_Data.handles(3));
    AMB_Data.handles(3)=0;
  end
end

% enough points
if np < AMB_Data.degree+1 
  if AMB_Data.handles(2)~=0  
    delete(AMB_Data.handles(2));
    AMB_Data.handles(2)=0;
  end
end
if np < 1 
  if AMB_Data.handles(1)~=0  
    delete(AMB_Data.handles(1));
    AMB_Data.handles(1)=0;
  end
end

% control polygon
if np>0
  if AMB_Data.handles(1)==0  
    AMB_Data.handles(1)=plot(AMB_Data.c(:,1),AMB_Data.c(:,2),'-ok');
  else
    set(AMB_Data.handles(1),'Xdata',AMB_Data.c(:,1));
    set(AMB_Data.handles(1),'Ydata',AMB_Data.c(:,2));
  end
end

% spline
if np > AMB_Data.degree
  if AMB_Data.movetype
    t=linspace(tau0,tau1,51);
  else
    t=linspace(tau0,tau1,301);
  end  
  b1=spline_evaluation(AMB_Data.c(:,1),AMB_Data.tau,t');
  b2=spline_evaluation(AMB_Data.c(:,2),AMB_Data.tau,t');
  if AMB_Data.handles(2)==0  
    AMB_Data.handles(2)=plot(b1,b2,'-b','linewidth',2);
  else
    set(AMB_Data.handles(2),'Xdata',b1);
    set(AMB_Data.handles(2),'Ydata',b2);
  end
end

if np < 1
  set(AMB_Data.uicontrols(2),'String','Mouse: left = insert point, right =  move knot ');
elseif np < AMB_Data.degree+1
  set(AMB_Data.uicontrols(2),'String','Mouse: left = insert, middle = delete point, right =  move point/knot ');
else
  set(AMB_Data.uicontrols(2),'String','Mouse: left = insert, right =  move point/knot, middle delete point ');
end

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

function p = spline_evaluation(c,xi,x) 
% SPLINE_EVALUATION: evaluation of a spline function
% p = spline_evaluation(c,xi,x) 
% values p(j) at x(j) of a spline function  
% with coefficients c and knots xi

% change to column vectors (if necessary)
c = c(:); xi = xi(:); x = x(:);

% degree
n = length(xi)-length(c)-1;

% add points and knots to avoid exceptional cases
c = [zeros(n+1,1); c; zeros(n+1,1)];
xi_max = (2*norm([xi; x],inf)+1)*ones(n+1,1);
xi = [-xi_max; xi; xi_max];

for j = 1:length(x)
   % interval: xi(kj) <= x(j) < xi(kj+1)
   kj = min(find(xi>x(j)))-1;
   % relevant control points
   pxj = c(kj-n:kj);
   % triangular scheme
   for m = 1:n 
      ind = [kj-n+m:kj];
      gamma = ((x(j) - xi(ind)) ./ (xi(ind+n+1-m) - xi(ind)));
      pxj = gamma.*pxj(2:end) + (1-gamma).*pxj(1:end-1);
   end 
   p(j) = pxj;
end
p = p(:); 

function [c,tau] = knot_insertion(c,tau,t) 
% KNOT_INSERTION: knot insertion for splines
% [c,tau] = knot_insertion(c,tau,t) 
% change control points c(k,:) and knots tau
% by adding knots t(j) in interior knot intervals with
% multiplicity at most equal to the degree (=> n>0)

% convert to column vectors (if necessary)
tau = tau(:); t = t(:);

% number of control points, dimension, and degree 
[m,dim] = size(c); n = length(tau) - m - 1;

% add knots successively 
for tj = t'
   % relevant indices
   ind = find(tau(2:end-n-1)<tj & tau(n+2:end-1)>tj)+1;
   if length(ind)>0
      % convex combinations
      gamma = repmat((tj-tau(ind))./(tau(ind+n)-tau(ind)),1,dim);
     c = [c(1:min(ind-1),:); 
          (1-gamma).*c(ind-1,:)+gamma.*c(ind,:); 
          c(max(ind):end,:)];
     tau = sort([tau; tj]);
   else 
      display(['knot ' num2str(tj) ' not admissible'])
   end
end

