procrustes and affine transform in Seg3D

[jessdtate]

This was suggested by Josh Blauer. The suggestion is to add a procrustes and/or affine transform. These algorithms are pretty straight forward, but require correspondence points. These could be chosen with the seed points widget. I would also be nice to be able to save and load the correspondence points (issue #46) and the transform. For the following algorithms, DATA and MODEL are 3xN matrices of the correspondence points. The output transforms are R (rotation),S(scaling), and T(translation) so that S_R_DATA+T(repeated for size(DATA)) = registered to MODEL. The transform to save would be [R*S T; 0 0 0 1]. The transformed label masked would have to be mapped onto a grid. Register point sets will probably have one of these algorithms implemented already and iterate to find the best match.

This is the matlab code for the constraint procrustes (rigid, no reflection):

med=mean(DATA,2); mem=mean(MODEL,2);
A=DATA-repmat(med,1,N);
B=MODEL-repmat(mem,1,N);

normA=sqrt(trace(A_A'));
normB=sqrt(trace(B_B'));
A=A/normA;
B=B/normB;

[U,Sig,V] = svd(B_A');
U(:,end)=U(:,end)_det(U_V');
R=U_V';

trsqrtAB = sum(diag(Sig));
S=trsqrtAB*normB/normA;

(There are some constraints on S and S can be made to be directional (scaling in the x y x directions)

This is the matlab code for affine transform:
P = zeros(3_N, 12);
q = zeros(3_N, 1);
ind = 1:3;
for n = 1:N
P(ind, :) = [DATA(1,n) 0 0 DATA(2,n) 0 0 DATA(3,n) 0 0 1 0 0; ....
0 DATA(1,n) 0 0 DATA(2,n) 0 0 DATA(3,n) 0 0 1 0;
0 0 DATA(1,n) 0 0 DATA(2,n) 0 0 DATA(3,n) 0 0 1 ];

q(ind) = MODEL(1:3, iclosest(n));
ind = ind + 3;
end;

theta=P\q;
R=[theta(1) theta(4) theta(7) ; theta(2) theta(5) theta(8); ...
theta(3) theta(6) theta(9)];
T=[theta(10); theta(11); theta(12)];
S=1;