function demo_spline
% graphical input of a spline 

global AMB_Data;
AMB_Data=[];
AMB_Data.degree=4;
AMB_Data.dimension=10;
AMB_Data.knots=linspace(0,1,AMB_Data.dimension+AMB_Data.degree+1)';
AMB_Data.c=[[5:2:23]'/28,ones(10,1)/2];
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.axis=[0,0];
AMB_Data.uicontrols=[];

% Initialize GUI
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',...
	  '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=move knots/coefficients ',...
	'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 dimension
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.25,.85,.2,.05],...
	'String','Dimension:',...
	'Style','text');

% Input for dimension
AMB_Data.uicontrols(2)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','right', ...
	'Position',[.4,.86,.06,.05],...
	'String',sprintf('%1d',AMB_Data.dimension),...
	'Style','edit',...
	'Callback',@dimension_update);

% 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;
myplot;


%%% Callback functions %%%

% Mouse pressed
function myDownfct(varargin)
global AMB_Data;

if strcmp(get(gcf,'SelectionType'),'normal') 
  p=get(AMB_Data.axis(1),'CurrentPoint'); % check for point
  p=p(1,1:2);  
  if min(p) >= 0 & max(p) <=1 
    % get 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=2; 
      AMB_Data.moveind=ind;
      AMB_Data.c(ind,2)=p(2);
      myplot;      
    end
  else  
    p=get(AMB_Data.axis(2),'CurrentPoint'); % check for knot
    p=p(1,1:2);  
    if p(1) >= 0 & p(1) <=1 & p(2)>-.1 & p(2) < .1 
      if ~isempty(AMB_Data.knots) 
	% nearest knot
	[m,ind]=min(abs(p(1)-AMB_Data.knots));
	if m < .05 % really near 
	  AMB_Data.movetype=1; 
	  AMB_Data.knots(ind)=p(1);
	  [AMB_Data.knots,sind]=sort(AMB_Data.knots);
	  AMB_Data.moveind=sind(ind);	  
	  myplot;      
	end      	
      end  
    end
  end
end

% Mouse released
function myUpfct(varargin)
global AMB_Data;
AMB_Data.movetype=0;
AMB_Data.moveind=0;
myplot;      

% Mouse movement
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;
  knots=[0;AMB_Data.knots;1];
  if p < knots(ind), p=knots(ind); end
  if p > knots(ind+2), p=knots(ind+2); end 
  AMB_Data.knots(ind)=p;
  
  grv=AMB_Data.knots(2:end-AMB_Data.degree);
  for k=2:AMB_Data.degree
    grv=grv+AMB_Data.knots(1+k:end-AMB_Data.degree+k-1);
  end
  grv=grv/AMB_Data.degree;
  AMB_Data.c(:,1)=grv;
  myplot;
elseif AMB_Data.movetype == 2
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,2);
  ind = AMB_Data.moveind;
  if p < 0, p=0; end
  if p > 1, p=1; end 
  AMB_Data.c(ind,2)=p;
  myplot;
end     

% New B-spline button 
function restart(varargin)
global AMB_Data;
knumber=AMB_Data.dimension+AMB_Data.degree+1;
AMB_Data.knots=linspace(0,1,knumber)';
AMB_Data.c=[[AMB_Data.degree+1:2:(2*knumber-3-AMB_Data.degree)]'/(2*knumber-2),...
		ones(AMB_Data.dimension,1)/2];
AMB_Data.movetype=0;
AMB_Data.moveind=0;
myplot;

% degree input
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

% dimension input
function dimension_update(varargin)
global AMB_Data;
dimension= str2double(get(AMB_Data.uicontrols(2),'String'));
if dimension < 1 | dimension > 25
  set(AMB_Data.uicontrols(2),'String',sprintf('%1d',AMB_Data.dimension));
else
  AMB_Data.dimension =dimension;
  restart;
end

%%% Internal Functions %%%

function myplot
global AMB_Data;
knots=AMB_Data.knots;

% 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;
% plot D_xi
xi0=knots(AMB_Data.degree+1);
xi1=knots(end-AMB_Data.degree);

if xi1 <= xi0 
  xiempty=1;
else
  xiempty=0;
  fill([xi0,xi0,xi1,xi1],[-.06,.06,.06,-.06],.5*[1,1,1]);  
end

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

% plot knots
n=size(AMB_Data.knots,1);
for k=1:length(knots) 
  plot(knots(k)*[1 1],[-.08,.08],'-k');
  if mult(k)>1
    text(knots(k)-.002,.12,sprintf('%d',mult(k)));
  end
end    
if AMB_Data.movetype  == 1
  plot(AMB_Data.knots(AMB_Data.moveind)*[1 1],[-.08,.08],'-r');
end  

% plot spline
subplot(AMB_Data.axis(1));
knots=AMB_Data.knots;
cla;
fill([-.1,1.1,1.1,-.1],[-.1,-.1,1.1,1.1],[1,1,1],'edgecolor','none');
hold on;
plot([0,1,1,0,0],[0,0,1,1,0],'-k');

if xiempty
  text(.5,.5,'Parameter interval empty',...
       'horizontalalign','center','fontweight','bold','fontsize',10);
  return;
end

if any(mult > AMB_Data.degree) | ...
      knots(AMB_Data.degree+1) >= knots(AMB_Data.degree+2) | ...
      knots(end - AMB_Data.degree-1) >= knots(end-AMB_Data.degree)
  text(.5,.5,sprintf('Knot sequence not %d-regular',AMB_Data.degree),...
       'horizontalalign','center','fontweight','bold','fontsize',10) ;
  return;
end

if n>1 
  if AMB_Data.movetype
    t=linspace(xi0,xi1,51);
  else
    t=linspace(xi0,xi1,301);
  end
  b=spline_evaluation(AMB_Data.c(:,2),AMB_Data.knots,t');
  plot(t,b,'-b','linewidth',2);
end

plot([0,1],1/3*[1,1],'-k');

% plot c_k
for k=1:size(AMB_Data.c,1) 
  plot(AMB_Data.c(k,1)*[1 1],[1/3,AMB_Data.c(k,2)],'-k','linewidth',2);
  plot(AMB_Data.c(k,1),AMB_Data.c(k,2),'ok');
end    

if AMB_Data.movetype  == 2
  plot(AMB_Data.c(AMB_Data.moveind,1)*[1 1],[1/3,AMB_Data.c(AMB_Data.moveind,2)],'-r','linewidth',2);
  plot(AMB_Data.c(AMB_Data.moveind,1),AMB_Data.c(AMB_Data.moveind,2),'or');
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(:); 

