function demo_rational_bezier_curve
% graphical input of a control polygon
% and weights for a rational Bezier curve

% Initialise data
global AMB_Data;
AMB_Data=[];
AMB_Data.plotpoints=300;
AMB_Data.movepoints=51;
AMB_Data.c=zeros(0,2);
AMB_Data.PW=zeros(0,2);
AMB_Data.w=[];
AMB_Data.handles=[0 0 0 0];
AMB_Data.movetype=0;
AMB_Data.moveind=0;
AMB_Data.uicontrols=[];
AMB_Data.axis=0;

% 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,.91,.9,.07],...
	'String','Rational Bezier Curve',...
	'Style','text');

% Text for degree
AMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','left', ...
	'ListboxTop',0, ...
	'Position',[.05,.84,.14,.05],...
	'String','Degree: -1', ...
	'Style','text');

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

% Text for click
uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'Fontweight','bold', ...
	'HorizontalAlignment','right', ...
	'Position',[.0,.01,.95,.05],...
	'String','Mouse: left = insert, middle = delete point, right = move point/weight ', ...
	'Style','text');


% Axis for curve
AMB_Data.axis(1)=subplot('position',[.05,.07,.9,.77]);
axis([-.01,1.01,-.01,1.01]);
axis manual;
axis off
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');

%%% Callback Functions %%%

% Reset button pressed
function restart(varargin) 
global AMB_Data;
AMB_Data.c=zeros(0,2);
AMB_Data.w=zeros(0,1);
AMB_Data.movetype=0;
AMB_Data.moveind=0;
for k=AMB_Data.handles
  if k,delete(k);end
end
AMB_Data.handles=[0,0,0,0];
set(AMB_Data.uicontrols(1),'String','Degree: -1');    

% Mouse button pressed
function myDownfct(varargin)
global AMB_Data;

if strcmp(get(gcf,'SelectionType'),'normal') % left mouse: new point 
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);  
  if min(p) >= 0 & max(p) <=1
    AMB_Data.c=[AMB_Data.c;p];
    AMB_Data.w=[AMB_Data.w;1];
    n=size(AMB_Data.c,1);
    AMB_Data.movetype=1;
    AMB_Data.moveind=n;
  end
  myplot;
end

if strcmp(get(gcf,'SelectionType'),'alt') % right mousse: change weight
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);  
  if min(p) >= 0 & max(p) <=1 & size(AMB_Data.c,1) >1 
    [m,ind]=min((p(1)-AMB_Data.PW(:,1)).^2+(p(2)-AMB_Data.PW(:,2)).^2);
    if m < .003 
      AMB_Data.movetype=4; % change weight ratio
      AMB_Data.moveind=ind;
      myplot;              
    else
      [m,ind]=min((p(1)-AMB_Data.c(:,1)).^2+(p(2)-AMB_Data.c(:,2)).^2);
      if m < .003 
	AMB_Data.movetype=1; % move point
	AMB_Data.moveind=ind;
	myplot;          
      end    
    end
  end
end

if strcmp(get(gcf,'SelectionType'),'extend') % middle mouse: delete point 
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);  
  if min(p) >= 0 & max(p) <=1 
    [m,ind]=min((p(1)-AMB_Data.c(:,1)).^2+(p(2)-AMB_Data.c(:,2)).^2);
    if m < .01 
      AMB_Data.c(ind,:)=[];
      AMB_Data.w(ind,:)=[];
      myplot;      
    end
  end
end

% Mouse button released
function myUpfct(varargin)
global AMB_Data;
AMB_Data.movetype=0;
AMB_Data.moveind=0;
if AMB_Data.handles(4)~=0
  set(AMB_Data.handles(4),'Xdata',NaN);
  set(AMB_Data.handles(4),'Ydata',NaN);  
  myplot;
end

% Mouse movement
function myMovefct(varargin)
global AMB_Data;
switch AMB_Data.movetype
 case 1 % move point
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);
  if min(p) >= 0 & max(p) <=1 
    AMB_Data.c(AMB_Data.moveind,:)=p;
    myplot;
  end 
 case 4
  ind=AMB_Data.moveind;
  P=AMB_Data.c;
  p=get(AMB_Data.axis(1),'CurrentPoint');
  p=p(1,1:2);
  tw=((p-P(ind,:))*diff(P(ind:ind+1,:))')/norm(diff(P(ind:ind+1,:)))^2;  
  if tw < .05, tw=.05;end
  if tw > .95, tw=.95;end
  AMB_Data.w(ind+1:end)=tw/(1-tw)*AMB_Data.w(ind)/AMB_Data.w(ind+1)*...
      AMB_Data.w(ind+1:end);
  % scale weights to sum n
  AMB_Data.w=AMB_Data.w*size(AMB_Data.w,1)/sum(AMB_Data.w);  
  myplot;
end     

%%% Internal Funcions %%%

function myplot
global AMB_Data;

if AMB_Data.movetype
  panz = AMB_Data.movepoints;
else
  panz = AMB_Data.plotpoints;
end  

% Update text for degree 
subplot(AMB_Data.axis(1));
n=size(AMB_Data.c,1);
set(AMB_Data.uicontrols(1),'String',sprintf('Degree: %1d',n-1));

% delete curve and control polygon 
if n < 2 
  if AMB_Data.handles(2)~=0  
    delete(AMB_Data.handles(2));
    AMB_Data.handles(2)=0;
  end
end
if n < 1 
  if AMB_Data.handles(1)~=0  
    delete(AMB_Data.handles(1));
    AMB_Data.handles(1)=0;
  end
end

if n>0
  % draw polygon
  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
if n>1  
  % draw curve 
  t=linspace(0,1,panz);
  [s,PW]=rational_bezier_curve(AMB_Data.c,AMB_Data.w,t');
  AMB_Data.PW=PW;
  if AMB_Data.handles(2)==0  
    AMB_Data.handles(2)=plot(s(:,1),s(:,2),'-b','Linewidth',2);  
  else
    set(AMB_Data.handles(2),'Xdata',s(:,1));
    set(AMB_Data.handles(2),'Ydata',s(:,2));
  end
  % draw weight points 
  if AMB_Data.handles(3)==0  
    AMB_Data.handles(3)= ...
	plot(AMB_Data.PW(:,1),AMB_Data.PW(:,2),'sm');    
  else    
    set(AMB_Data.handles(3),'Xdata',AMB_Data.PW(:,1));
    set(AMB_Data.handles(3),'Ydata',AMB_Data.PW(:,2));
  end  
end

% update marked weight point 
if AMB_Data.movetype==4
  if AMB_Data.handles(4)==0  
    AMB_Data.handles(4)=...
	plot(AMB_Data.PW(AMB_Data.moveind,1),AMB_Data.PW(AMB_Data.moveind,2),...
	     'om','markersize',20);          
  else
    set(AMB_Data.handles(4),'Xdata',AMB_Data.PW(AMB_Data.moveind,1));
    set(AMB_Data.handles(4),'Ydata',AMB_Data.PW(AMB_Data.moveind,2));
  end
end

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

function [r,d] = rational_bezier_curve(c,w,t) 
% RATIONAL_BEZIER_CURVE: evaluation of a rational Bezier curve 
% [r,d] = rational_bezier_curve(c,w,t) 
% points r(j,:) at t(j) and weight points d(i,:) 
% for a rational Bezier curve 
% with control points c(k,:) and weights w

% convert to column vector (if necessary)
w = w(:);

% degree and dimension
[np1,dim] = size(c);

% homogeneous control points
W = repmat(w,1,dim);
cw = [W.*c w];

% evaluation
r = casteljau(cw,t);
r= r(:,1:end-1)./repmat(r(:,end),1,dim);

% weight points
K = [1:np1-1];
d = (W(K,:).*c(K,:)+W(K+1,:).*c(K+1,:))./(W(K,:)+W(K+1,:));

function p = casteljau(c,t)
% CASTELJAU: de Casteljau's algorithm for Bezier curves 
% p = casteljau(c,t)
% points p(j,:) on a Bezier curve with control points c(k,:)
% for parameters t(j) 

% convert to column vector (if necessary)
t = t(:); 
 
[np1,dim] = size(c);
% de Casteljau's recursion 
if np1 == 1
   p = (1+0*t)*c;
else
   p0 = casteljau(c(1:end-1,:),t); 
   p1 = casteljau(c(2:end,:),t); 
   t = repmat(t,1,dim); 
   p = (1-t).*p0 + t.*p1;
end

