function [handle, rhandle] = expVisualiseRootProj(channels, skel, padding)

% Version 0.100 (Unsupported, unreleased) 
%
% Code provided by Graham Taylor and Geoff Hinton 
%
% For more information, see:
%     http://www.cs.toronto.edu/~gwtaylor/publications/nips2006mhmublv
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, expressed or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.
%
% Based on skelVisualise.m version 1.1
% Copyright (c) 2006 Neil D. Lawrence
%
% We support two types of skeletons:
%  1) Those built from the CMU database (acclaim)
%     http://mocap.cs.cmu.edu/
%  2) Those built from data from Eugene Hsu (mit)
%     http://people.csail.mit.edu/ehsu/work/sig05stf/
% EXPVISUALISE For updating a skel representation of 3-D data.
%
% Usage: handle = expVisualise(channels, skel, padding)


%addpath ../Toolboxes/rotationMatrix


 

if nargin<3
    padding = 0;
end
channels = [channels zeros(1, padding)];
vals = exp2xyz(skel, channels);
connect = skelConnectionMatrix(skel);

indices = find(connect);
[I, J] = ind2sub(size(connect), indices);

handle(1) = plot3(vals(:, 1), vals(:, 3), vals(:, 2), '.');

%this is certainly not needed for MIT data
if any(strcmp(skel.type,{'acclaim','bvh','cmubvh'}))
  axis ij % make sure the left is on the left. Added in MOCAP0p134
end

set(handle(1), 'markersize', 20);
hold on
grid on
for i = 1:length(indices)
    handle(i+1) = line([vals(I(i), 1) vals(J(i), 1)], ...
        [vals(I(i), 3) vals(J(i), 3)], ...
        [vals(I(i), 2) vals(J(i), 2)]);
    set(handle(i+1), 'linewidth', 2);
end


d0 = [0 0 1]';
l0 = [1 0 0]';
v0 = [0 1 0]';


 expRot = channels(skel.tree(1).expmapInd);
 %expRot = channels([5 6 4]);
 Rot = RotationMatrix(expRot,'exponentialMap');
 
 R = Rot.GetRotationMatrix();
 
 D = R' * d0;
 L = R' * l0;
 A = R' * v0;  % we don't use this!
 
 D = D*20;
 L = L*20;
 A = A*20;
 
 
  %roll
  phi = acos(dot(D,-v0)/norm(D)); %this corresponds to V in the formulas
  
  %pitch
  psi = acos(dot(L,-v0)/norm(L)); %this corresponds to U in the formulas
  
  %yaw
  top = find(abs(d0)==1);
  bot = find(abs(l0)==1);
  
  
  theta = atan2(L(top),L(bot));   %% be careful: choose the components of the u based on the index d0 and l0 point
  
  
  Dxz = [D(3) D(1) 0];
  Lxz = [L(3) L(1) 0];
  
 rhandle(1) = quiver3(vals(1, 1), vals(1, 3), vals(1, 2),Dxz(1),Dxz(2),Dxz(3),'linewidth', 4, 'color','[0.5,0.0,0.0]'); %u
 rhandle(2) = quiver3(vals(1, 1), vals(1, 3), vals(1, 2),Lxz(1),Lxz(2),Lxz(3),'linewidth', 4,'color','[0.0,0.5,0]'); %v
 %rhandle(3) = quiver3(vals(1, 1), vals(1, 3), vals(1, 2),w(1),w(2),w(3),'linewidth', 4,'color','[0.0,0.0,0.5]'); %w
 
axis equal
xlabel('x')
ylabel('z')
zlabel('y')
axis on

end
