% Hands.m
% 
% Using Ramsay's "handwrit.m" example as inspiration for registration of 
% a dataset of 40 hands, each consisting og 200 landmarks. 
% Plots landmarks, fd-curves, velocity and acceleration. 
% B-splines and Fourier series are applied.
%
% Søren G. Erbou
% 29/3-05
% Phd-course 02915

close all
clear all
clc

addpath ('E:\SGE\PhD\02915\FDA\fdaM')

WindowPos=get(0,'screensize');
WindowPos=[50,50,0.9*WindowPos(3),0.85*WindowPos(4)];
FontSize=14;

load hands
[Y,X,Z]=size(hands); % [200,40,2]
fdarray = hands;

% Parameterization of landmarks
fdaparam   = (1:Y)';

% Smoothing parameter
lambda   = 1e0;

Basis=1; % 1=B-spline, 2=Fourier

switch Basis
    case 1
    % Bspline setup
    par(1)=101; % No of basis functions
    par(2)=6;   % Order of B-splines
    fdabasisobj = create_bspline_basis([1, Y], par(1), par(2));
    Basis='B-spline';
    case 2
    % Fourier setup
    par(1)=64; % No of basis functions
    fdabasisobj = create_fourier_basis([1, Y], par(1));
    Basis='Fourier';
end
npar=length(par);

%  set up the functional data structure without smoothing

fdafd = data2fd(fdarray, fdaparam, fdabasisobj);
fdafd_fdnames{1} = 'Landmarks';
fdafd_fdnames{2} = 'Replications';
fdafd_fdnames{3} = 'Units';
fdafd = putnames(fdafd, fdafd_fdnames);

%  compute values of curves and the values of the mean curve

fdamat     = eval_fd(fdafd, fdaparam);
fdameanvec = squeeze(eval_fd(mean(fdafd), fdaparam));

%  plot individual curves, including both sampling points and fit
%  also plot the mean curve in the background

F1=figure;
set(F1,'Position',WindowPos);
for i = 1:40
    plot(fdarray(:,i,1),  fdarray(:,i,2), 'go', ...
        fdamat(:,i,1),   fdamat(:,i,2), '-', ...
        fdameanvec(:,1), fdameanvec(:,2), 'r--');
%     hold on
    axis equal
    axis([0.2 1.25 0.1 1.0]);
    legend('Landmarks','FD curve','Mean shape');
    xlabel('X','FontSize',FontSize)
    ylabel('Y','FontSize',FontSize)
    Title=['Hand no. ', num2str(i),', ',Basis,', ',num2str(par(1))];
    j=2;
    while j<=npar
        Title=[Title,', ',num2str(par(j))];
        j=j+1;
    end
    title(Title,'FontSize',FontSize);
    pause(0.1);
end
pause
%  apply a small amount of smoothing to get better acceleration
%  curves

%  set up the functional parameter object to define amount of
%  smoothing

Lfdobj   = int2Lfd(2);
fdParobj = fdPar(fdafd, Lfdobj, lambda);

%  smooth the functional data object

 fdasmthfd = smooth_fd(fdafd, fdParobj);

%  plot the velocity records
F2=figure;
set(F2,'Position',WindowPos);
for i=1:40
    veloc = squeeze(eval_fd(fdaparam, fdasmthfd(i,:), int2Lfd(1)));
    veloc_rough = squeeze(eval_fd(fdaparam, fdafd(i,:), int2Lfd(1)));
    subplot(2,1,1)
    plot(fdaparam, veloc_rough(:,1), '-', fdaparam, veloc_rough(:,2), '-', ...
        [1,200],[0,0], 'r:');
    %     hold on
    axis([1 200 -0.025 0.025]);
    xlabel('Landmarks','FontSize',FontSize)
    ylabel('Units/Landmarks','FontSize',FontSize)
    title([num2str(i),' - 1st deriv. - no smoothing'],'FontSize',FontSize)
    legend('X','Y','Location','NorthWest')
    subplot(2,1,2)
    plot(fdaparam, veloc(:,1), '-', fdaparam, veloc(:,2), '-', ...
        [1,200],[0,0], 'r:');
    %     hold on
    xlabel('Landmarks','FontSize',FontSize)
    ylabel('Units/Landmarks','FontSize',FontSize)
    title([,num2str(i),' - 1st deriv. - smoothed, \lambda = ',num2str(lambda,'%2.0e')],'FontSize',FontSize)
    legend('X', 'Y','Location','NorthWest')
    axis([1 200 -0.025 0.025]);
    pause(0.2)
end
pause

F3=figure;
set(F3,'Position',WindowPos);
for i=1:40
    accel = squeeze(eval_fd(fdaparam, fdasmthfd(i,:), int2Lfd(2)));
    accel_rough = squeeze(eval_fd(fdaparam, fdafd(i,:), int2Lfd(2)));
    subplot(2,1,1)
    plot(fdaparam, accel_rough(:,1), '-', fdaparam, accel_rough(:,2), '-', ...
        [1,200],[0,0], 'r:');
    %     hold on
    axis([1 200 -0.02 0.02]);
    xlabel('Landmarks','FontSize',FontSize)
    ylabel('Units/Landmarks/Landmarks','FontSize',FontSize)
    title([num2str(i),' - 2nd deriv. - no smoothing'],'FontSize',FontSize)
    legend('X','Y','Location','NorthWest')
    subplot(2,1,2)
    plot(fdaparam, accel(:,1), '-', fdaparam, accel(:,2), '-', ...
        [1,200],[0,0], 'r:');
    %     hold on
    xlabel('Landmarks','FontSize',FontSize)
    ylabel('Units/Landmarks/Landmarks','FontSize',FontSize)
    title([,num2str(i),' - 2nd deriv. - smoothed, \lambda = ',num2str(lambda,'%2.0e')],'FontSize',FontSize)
    legend('X', 'Y','Location','NorthWest')
    axis([1 200 -0.02 0.02]);
    pause(0.2)
end
pause
%  compute the acceleration magnitudes

D2fdamat = eval_fd(fdasmthfd, fdaparam, int2Lfd(2));
D2mag = sqrt(D2fdamat(:,:,1).^2 + D2fdamat(:,:,2).^2);

%  plot the acceleration magnitudes

F4=figure;
set(F4,'Position',WindowPos);
plot(fdaparam, D2mag)
xlabel('Landmarks','FontSize',FontSize)
ylabel('Units/Landmarks/Landmarks','FontSize',FontSize)
title('Acceleration magnitude','FontSize',FontSize);
axis([1 200 0 0.02]);