Source Code for Module geometry.procrustes

 1  from . import contract, np 
 2   
 3  from contracts import check 
 4   
 5   
 6  # TODO: write tests 
 7  @contract(X='array[KxN],K>=2,K<N', Y='array[KxN]', 
 8            returns='array[KxK],orthogonal') 
 9 -def best_orthogonal_transform(X, Y): 
10      ''' Finds the best orthogonal transform R  between X and Y, 
11          such that R X ~= Y. ''' 
12      YX = np.dot(Y, X.T) 
13      check('array[KxK]', YX) 
14      U, _, V = np.linalg.svd(YX) 
15      best = np.dot(U, V) 
16      return best 
17   
18   
19  @contract(M='array[NxN]', returns='array[NxN],orthogonal') 
20 -def closest_orthogonal_matrix(M): 
21      ''' Finds the closest orthogonal matrix to M. ''' 
22      U, _, V = np.linalg.svd(M) 
23      R = np.dot(U, V) 
24      return R 
25   
26   
27  # TODO: write tests 
28  @contract(X='array[KxN],K>=2,K<N', Y='array[KxN]', 
29            returns='tuple( (array[KxK],orthogonal), array[Kx1])') 
30 -def best_similarity_transform(X, Y): 
31      ''' Finds the best transform (R,t)  between X and Y, 
32          such that R X + t ~= Y. ''' 
33      K = X.shape[0] 
34      Xm = X.mean(axis=1).reshape(K, 1) 
35      Ym = Y.mean(axis=1).reshape(K, 1) 
36      X = X - Xm 
37      Y = Y - Ym 
38  #    assert_allclose(X.mean(axis=1), 0, atol=1e-8) 
39  #    assert_allclose(Y.mean(axis=1), 0, atol=1e-8) 
40      YX = np.dot(Y, X.T) 
41      check('array[KxK]', YX) 
42      U, _, V = np.linalg.svd(YX) 
43      R = np.dot(U, V) 
44      t = Ym - np.dot(R, Xm) 
45      return R, t 
46