function demo_bezier_patch
% morphing of bezier patches

% Initialize data
global AMB_Data;
AMB_Data=[];
AMB_Data.c{1}=cat(3,...
		 [0 0 0 0;1 1 1 1;2 2 2 2;3 3 3 3],...
		 [0 0 0 0;0 0 0 0;3 3 3 3;3 3 3 3]',...
		 [0 1 1 0;0 1 1 0;0 1 1 0;0 1 1 0]*4/3);

AMB_Data.c{2}=cat(3,...
		 [0 0 0 0;1 1 1 1;2 2 2 2;3 3 3 3],...
		 [0 1 1 0;0 1 1 0;3 2 2 3;3 2 2 3]',...
		 [0 2 2 0;0 1 1 0;0 1 1 0;0 2 2 0]*2/3);


AMB_Data.c{3}=cat(3,...
		 [0 0 0 0;1 1 1 1;2 2 2 2;3 3 3 3],...
		 [0 0 0 0;1 1 1 1;2 2 2 2;3 3 3 3]',...
		 [0 1 1 0;0 1 1 0;0 -1 -1 0;0 -1 -1 0]*4/3);

AMB_Data.c{4}=cat(3,...
		 [0 0 0 0;0 -1 -1 0;0 1 1 0;0 0 0 0],...
		 [0 0 0 0;0 -1 -1 0;0 1 1 0;0 0 0 0]',...
		 [0 0 0 0;0 -1 -1 0;0 -1 -1 0;0 0 0 0]);

% Initialize GUI
clf reset;
% Text for title
uicontrol('Units','normalized', ...
  	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','center', ...
        'FontWeight','bold',...
        'FontSize',20,...
	'Position',[.05,.91,.9,.07],...
	'String','Bezier Patch',...
	'Style','text');

% Popup for source
AMB_Data.uicontrols(1)=uicontrol('Units','normalized', ...
	  'BackgroundColor',get(gcf,'Color'), ...
	  'HorizontalAlignment','left', ...
	  'Value',1,...
	  'Position',[.3,.85,.2,.05],...
          'Callback',@plotsource,...
	  'String','tunnel|bridge|slide|drop',...
	  'Style','popup');

uicontrol('Units','normalized', ...
  	'BackgroundColor',get(gcf,'Color'), ...
	'HorizontalAlignment','center', ...
        'FontSize',10,...
	'Position',[.60,.84,.05,.05],...
	'String','to',...
	'Style','text');


% Popup for destination 
AMB_Data.uicontrols(2)=uicontrol('Units','normalized', ...
	  'BackgroundColor',get(gcf,'Color'), ...
	  'HorizontalAlignment','left', ...
	  'Value',2,...
	  'Position',[.75,.85,.2,.05],...
	  'String','tunnel|bridge|slide|drop',...
	  'Style','popup');

% Button for morphing
AMB_Data.uicontrols(3)=uicontrol('Units','normalized', ...
	'BackgroundColor',get(gcf,'Color'), ...
	'String','Morph', ...
	'Enable','On', ...
	'Position',[.05,.85,.10,.05], ...
	'Callback',@morph,...
	'Style','pushbutton');

% Axis 
AMB_Data.axis(1)=subplot('position',[.05,.07,.9,.7]);
axis equal;
axis off;
cm=winter(1000);
colormap(cm(200:800,:))
view(51,26);
rotate3d on;
myplot(AMB_Data.c{1});

%%% Callback Functions %%%
function plotsource(varargin) 
global AMB_Data;
myplot(AMB_Data.c{get(AMB_Data.uicontrols(1),'value')});

function morph(varargin) 
global AMB_Data;

ind1=get(AMB_Data.uicontrols(1),'Value');
ind2=get(AMB_Data.uicontrols(2),'Value');
if ind1==ind2,return,end;
set(AMB_Data.uicontrols(3),'Enable','off');

for t=linspace(0,1,50);
  c=(1-t)*AMB_Data.c{ind1}+t*AMB_Data.c{ind2};
  myplot(c);
  pause(.1)
end
set(AMB_Data.uicontrols(1),'Value',ind2);
set(AMB_Data.uicontrols(2),'Value',mod(ind2,4)+1);
set(AMB_Data.uicontrols(3),'Enable','on');

%%% Internal Functions %%%
function myplot(c)

cla;
axis equal;
axis off;
cm=winter(1000);
colormap(cm(200:800,:))
view(51,26);
hold on;
mesh(c(:,:,1),c(:,:,2),c(:,:,3),...
     'facecolor','none','edgecolor','black','linewidth',1,...
     'markersize',10,'marker','o','markerfacecolor','black');
s=bezier_patch(c,linspace(0,1,51),linspace(0,1,51));
surf(s(:,:,1),s(:,:,2),s(:,:,3),'facecolor','interp','edgecolor','none','facealpha',.5);
plot3(s(1,:,1),s(1,:,2),s(1,:,3),'-k','linewidth',1.5);
plot3(s(end,:,1),s(end,:,2),s(end,:,3),'-k','linewidth',1.5);
plot3(s(:,1,1),s(:,1,2),s(:,1,3),'-k','linewidth',1.5);
plot3(s(:,end,1),s(:,end,2),s(:,end,3),'-k','linewidth',1.5);

plot3(s(26,:,1),s(26,:,2),s(26,:,3),'-k','linewidth',.5);
plot3(s(:,26,1),s(:,26,2),s(:,26,3),'-k','linewidth',.5);
axis image


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

function p = bezier_patch(c,t1,t2)
% BEZIER_PATCH: evaluation of a Bezier patch  
% p = bezier_patch(c,t1,t2)
% points p(j1,j2,:) of a Bezier patch with control points c(k1,k2,:)
% for parameters t1(j1) and t2(j2) 

[n1,n2,dim] = size(c);
for i=1:dim
   b1 = bernstein(n1-1,t1); b2 = bernstein(n2-1,t2); 
   p(:,:,i) = b1*c(:,:,i)*b2';
end

function b = bernstein(n,t)
% BERNSTEIN: Bernstein polynomials
% b = bernstein(n,t)
% values b(j,k) of b^n_{k-1} at t(j)

% change to column vector (if necessary)
t=t(:);

% initialization for degree 0 
b = zeros(length(t),2*n+1); 
b(:,n+1) = 1;

% recursion for Bernstein polynomials
for k=1:n
   tk = repmat(t,1,2*n+1-k);
   b = (1-tk).*b(:,2:end) + tk.*b(:,1:end-1);
end


