function [cx,cy,E,L]=bookstein(X,Y,beta_k);
% [cx,cy,E,L]=bookstein(X,Y,beta_k);
%
% Bookstein PAMI89

N=size(X,1);
Nb=size(Y,1);

if N~=Nb
	error('number of landmarks must be equal')
end

% compute distances between left points
r2=dist2_C(X,X);

% if min(min(r2+eye(N,N)))<.00000001	keyboard;	end

K=r2.*log(r2+eye(N,N)); % add identity matrix to make K zero on the diagonal
P=[ones(N,1) X];
L=[K  P
   P' zeros(3,3)];
V=[Y' zeros(2,3)];
if nargin>2
	% regularization
	L(1:N,1:N)=L(1:N,1:N)+beta_k*eye(N,N);
end

% if abs(rcond(L))<1e-15	keyboard; end

invL=inv(L);

c=invL*V';
cx=c(:,1);
cy=c(:,2);

if nargout>2
	% compute bending energy (w/o regularization)
	Q=c(1:N,:)'*K*c(1:N,:);
	E=mean(diag(Q));
end