function demo_spline_bivariate
% graphical movement of knots and
% input of a function for coefficients 
% for a bivariate spline

global AMB_Data;
AMB_Data=[];
AMB_Data.degree1=2;
AMB_Data.degree2=2;
AMB_Data.dimension1=10;
AMB_Data.dimension2=10;
AMB_Data.xi1=linspace(0,1,AMB_Data.dimension1+AMB_Data.degree1+1)';
AMB_Data.xi2=linspace(0,1,AMB_Data.dimension2+AMB_Data.degree2+1)';
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];
AMB_Data.Fun = @(x,y) .5*cos(2*pi*x)+sin(pi*y);  
AMB_Data.Funstring = '@(x,y) .5*cos(2*pi*x)+sin(pi*y)';  
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 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 = move knots ',...
	'Style','text');

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

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

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

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

% Text for dimensions
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.3,.85,.2,.05],...
	'String','Dimensions:',...
	'Style','text');

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

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

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

% Text for Function
uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.05,.80,.2,.05],...
	'String','Function:',...
	'Style','text');


% Input for Function
AMB_Data.uicontrols(7)=uicontrol('Units','normalized', ...
        'BackgroundColor',get(gcf,'Color'), ...
        'HorizontalAlignment','left', ...
	'Position',[.18,.80,.77,.05],...
	'String',sprintf('%s',AMB_Data.Funstring),...
	'Style','edit',...
	'Callback',@function_update);

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

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

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

%%% Callback Functions %%%

function myDownfct(varargin)
global AMB_Data;

if strcmp(get(gcf,'SelectionType'),'normal') % left mouse 
  p=get(AMB_Data.axis(2),'CurrentPoint'); % horizontal knot 
  p=p(1,1:2);  
  if p(1) >= 0 & p(1) <=1 & p(2)>-.1 & p(2) < .1 
    if ~isempty(AMB_Data.xi1) 
      % nearest knot
      [m,ind]=min(abs(p(1)-AMB_Data.xi1));
      if m < .05 
	AMB_Data.movetype=1; 
	AMB_Data.xi1(ind)=p(1);
	[AMB_Data.xi1,sind]=sort(AMB_Data.xi1);
	AMB_Data.moveind=sind(ind);	  
	myplot;      
      end      	
    end 
  else 
    p=get(AMB_Data.axis(3),'CurrentPoint'); % vertical knot
    p=p(1,1:2);  
    if p(2) >= 0 & p(2) <=1 & p(1)>-.1 & p(1) < .1 
      if ~isempty(AMB_Data.xi2) 
	[m,ind]=min(abs(p(2)-AMB_Data.xi2));
	if m < .05 
	  AMB_Data.movetype=2; 
	  AMB_Data.xi2(ind)=p(2);
	  [AMB_Data.xi2,sind]=sort(AMB_Data.xi2);
	  AMB_Data.moveind=sind(ind);	  
	  myplot;      
	end      	
      end 
    end
  end
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 = move knots ');
else
  AMB_Data.rotate=1;
  rotate3d on;
  set(AMB_Data.uicontrols(1),'String','disable rotation to change knots ');
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;
  xi=[0;AMB_Data.xi1;1];
  if p < xi(ind), p=xi(ind); end
  if p > xi(ind+2), p=xi(ind+2); end 
  AMB_Data.xi1(ind)=p;
  myplot;  
elseif  AMB_Data.movetype==2
  p=get(AMB_Data.axis(3),'CurrentPoint');
  p=p(1,2);
  ind = AMB_Data.moveind;
  xi=[0;AMB_Data.xi2;1];
  if p < xi(ind), p=xi(ind); end
  if p > xi(ind+2), p=xi(ind+2); end 
  AMB_Data.xi2(ind)=p;
  myplot;
end     

function restart(varargin)
global AMB_Data;
AMB_Data.xi1=linspace(0,1,AMB_Data.dimension1+AMB_Data.degree1+1)';
AMB_Data.xi2=linspace(0,1,AMB_Data.dimension2+AMB_Data.degree2+1)';
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.rotate=0;
rotate3d off;
set(AMB_Data.uicontrols(1),'String','Mouse: left = move knots ');
set(AMB_Data.uicontrols(6),'Value',0);
myplot;




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

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

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

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

% new function input
function function_update(varargin)
global AMB_Data;
AMB_Data.Funstring = get(AMB_Data.uicontrols(7),'String');
AMB_Data.Fun= eval(AMB_Data.Funstring);
restart;

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

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

% plot D_xi1
xi1=AMB_Data.xi1;
ximin=xi1(AMB_Data.degree1+1);
ximax=xi1(end-AMB_Data.degree1);

if ximax <= ximin 
  xiempty1=1;
else
  xiempty1=0;
  fill([ximin,ximin,ximax,ximax],[-.06,.06,.06,-.06],.5*[1,1,1]);
end

% calculate multiplicities
[xi,mult]=multiplicities(xi1);

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

% check regularity
irreg1 = irregular(xi1,mult,AMB_Data.degree1);

% set evaluatipn values
if AMB_Data.movetype
  x1=linspace(ximin,ximax,51)';
else
  x1=linspace(ximin,ximax,301)';
end

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

% plot D_xi2
xi2=AMB_Data.xi2;
ximin=xi2(AMB_Data.degree2+1);
ximax=xi2(end-AMB_Data.degree2);

if ximax <= ximin 
  xiempty2=1;
else
  xiempty2=0;
  fill([-.06,.06,.06,-.06],[ximin,ximin,ximax,ximax],.5*[1,1,1]);
end

% calculate multiplicities
[xi,mult]=multiplicities(xi2);

% plot knots xi2
for k=1:length(xi) 
  plot([-.08,.08],xi(k)*[1 1],'-k');
  if mult(k)>1
    text(.12,xi(k)-.002,sprintf('%d',mult(k)));
  end
end    

if AMB_Data.movetype ==2 
  plot([-.08,.08],AMB_Data.xi2(AMB_Data.moveind)*[1 1],'-r');
end  

% check regularity
irreg2 = irregular(xi2,mult,AMB_Data.degree2);

if AMB_Data.movetype
  x2=linspace(ximin,ximax,51)';
else
  x2=linspace(ximin,ximax,301)';
end

% spline
subplot(AMB_Data.axis(1));
cla;
axis off 

% check regularity
hold on;
if xiempty1 | xiempty2
  text(.5,.5,'Parameter interval empty',...
       'horizontalalign','center','fontweight','bold','fontsize',10);
  return;
end

if irreg1 | irreg2 
  text(.5,.5,'Knot sequence not regular',...
       'horizontalalign','center','fontweight','bold','fontsize',10) ;
  return;
end

% evaluate function at greville points
grv1=AMB_Data.xi1(2:end-AMB_Data.degree1);
for k=2:AMB_Data.degree1
  grv1=grv1+AMB_Data.xi1(1+k:end-AMB_Data.degree1+k-1);
end
grv1=grv1/AMB_Data.degree1;

grv2=AMB_Data.xi2(2:end-AMB_Data.degree2);
for k=2:AMB_Data.degree2
  grv2=grv2+AMB_Data.xi2(1+k:end-AMB_Data.degree2+k-1);
end
grv2=grv2/AMB_Data.degree2;
[c1,c2]=ndgrid(grv1,grv2);
c3=AMB_Data.Fun(c1,c2);

% plot spline
b=spline_evaluation_bivariate(c3,xi1,xi2,x1,x2);
surf(x1,x2,b','edgecolor','none','facecolor','interp','facealpha',.5);
hold on; 
mesh(c1,c2,0*c3,'facecolor','none','edgecolor','k');
plot3([1;1]*c1(:)',[1;1]*c2(:)',[0;1]*c3(:)','-k');
plot3(c1(:),c2(:),c3(:),'ok','Markersize',4);
axis image

function [xi,mult]=multiplicities(xi)
dxi=diff(xi);
ind=find(dxi<1000*eps);
mult=ones(size(xi));
for k=fliplr(ind')
  mult(k)=mult(k)+mult(k+1);
  mult(k+1)=[];
  xi(k+1)=[];
end

function irreg=irregular(xi,mult,deg)
if any(mult > deg) | xi(deg+1) >= xi(deg+2) | xi(end - deg-1) >= xi(end-deg)
  irreg=1;
else
  irreg=0;
end

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

function p = spline_evaluation_bivariate(c,xi1,xi2,x1,x2) 
% SPLINE_EVALUATION_BIVARIATE: evaluation of a bivariate spline function
% p = spline_evaluation_bivariate(c,xi1,xi2,x1,x2) 
% values p(j1,j2) of a spline with coefficients c(k1,k2)
% and knots xi1 and xi2 at (x1(j1),x2(j2)) 

% first parameter direction
c = spline_curve(c,xi1(:),x1(:));
% second parameter direction
p = spline_curve(c',xi2(:),x2(:))';


function p = spline_curve(c,tau,t) 
% SPLINE_CURVE: evaluation of a spline curve
% p = spline_curve(c,tau,t) 
% points p(j) for parameters t(j) of a spline curve 
% with control points c(k,:) and knots tau

% degree and dimension
[lc,dim] = size(c);
n = length(tau)-lc-1;

% add points and knots to avoid exceptional cases
c = [zeros(n+1,dim); c; zeros(n+1,dim)];
tau_max = repmat(2*norm([tau; t],inf),n+1,1);
tau = [-tau_max; tau; tau_max];


for j = 1:length(t)
  % interval: tau(kj) <= t(j) < tau(kj+1)
   kj = find(tau(1:end-1) <= t(j) & ... 
       tau(2:end) > t(j));
   % relevant control points
   ptj = c(kj-n:kj,:);
   % triangular scheme
   for m = 1:n 
      ind = [kj-n+m:kj];
      gamma = repmat(((t(j) - tau(ind)) ./ ... 
         (tau(ind+n+1-m) - tau(ind))),1,dim);
      ptj = gamma .* ptj(2:end,:) + ... 
         (1-gamma) .* ptj(1:end-1,:);
   end 
   p(j,:) = ptj;
end 






