function [C2, E] = thinplate(S, Q)

% thinplate(S, Q)
%
% Calcula um par de thin-plate splines que representam a
% transformação S -> Q.
%
% Recebe como argumentos um conjunto de marcos 
%
%     |s1_x s1_y|
% S = |s2_x s2_y|
%     |...  ... |
%     |sn_x sn_y|
%
% que representam a forma de entrada e um conjunto de pontos de
% controle
%
%     |q~1_x q~1_y|
% Q = |q~2_x q~2_y|
%     | ...   ... |
%     |q~n_x q~n_y|
%
% que representam a transformação de S por um par de thin-plate
% splines ou uma forma com a qual queremos comparar S.
%
% Devolve em C2 a matriz 
%
% C2 = |wx_1 wx_2 ... wx_n a_x b_xx b_xy|'
%      |wy_1 wy_2 ... wy_n a_y b_yx b_yy|
%
% contendo os parâmetros das thin-plate splines que representam a
% transformação S -> Q. 
%
% Em E é devolvido o valor da energia de dobramento correspondente.
%
% Referência:
%
% R. M. Cesar Jr. e L. F. Costa. "Shape Analysis and
% Classification". CRC Press, 2001, seção 4.9.6.

% Número de pontos de controle
n = size(Q, 1);

% Matriz T
for i=1:n
  for j=1:n
    T(i, j) = base(dist([S(i, 1), S(i, 2)], [S(j, 1), S(j, 2)]));
  end
end

% Matriz Se
for i=1:n
  Se(i, 1) = 1;
  Se(i, 2) = S(i, 1);
  Se(i, 3) = S(i, 2);
end

% Matriz M
M = [T, Se;
     Se', zeros(3)];

% Matriz H2
for i=1:n
  H2(i, 1) = Q(i, 1);
  H2(i, 2) = Q(i, 2);
end

for i=n+1:n+3
  H2(i, 1) = 0;
  H2(i, 2) = 0;
end

% Cálculo de C2
C2 = inv(M) * H2;

% Matriz W
for i=1:n
  W(i, 1) = C2(i, 1);
  W(i, 2) = C2(i, 2);
end

% Cálculo da energia de dobramento
E = trace(W' * T * W);



