Source Code for Module geometry.manifolds.special_orthogonal_group

 1  from . import DifferentiableManifold, MatrixLieGroup, np, S2, so, contract 
 2  from .. import (assert_allclose, rot2d, random_rotation, 
 3      axis_angle_from_rotation, rotation_from_axis_angle) 
 4  from contracts import check 
 5   
 6   
 7 -class SO_group(MatrixLieGroup): 
 8      '''  
 9          This is the Special Orthogonal group SO(n) describing rotations 
10          of Euclidean space; implemented for n=2,3. 
11           
12          TODO: do SO2 and SO3 separately 
13      ''' 
14   
15      @contract(N='int,(2|3)') 
16 -    def __init__(self, N): 
17          algebra = so[N] 
18          dimension = {2: 1, 3: 3}[N] 
19          MatrixLieGroup.__init__(self, n=N, algebra=algebra, 
20                                  dimension=dimension) 
21          DifferentiableManifold.embedding(self, algebra, 
22                                            self.algebra_from_group, 
23                                      self.group_from_algebra, 
24                                      itype='lie') 
25   
26 -    def __repr__(self): 
27          return 'SO%s' % (self.n) 
28   
29 -    def belongs(self, x): 
30          # TODO: make this much more efficient 
31          check('array[NxN],orthogonal', x, N=self.n) 
32          det = np.linalg.det(x) 
33          assert_allclose(det, 1, err_msg='I expect the determinant to be +1.') 
34   
35      @contract(returns='belongs') 
36 -    def sample_uniform(self): 
37          if self.n == 2: 
38              theta = np.random.rand() * np.pi 
39              return rot2d(theta) 
40          elif self.n == 3: 
41              return random_rotation() 
42          else: 
43              assert False, 'Not implemented for n>=4.' 
44   
45      @contract(x='belongs') 
46 -    def friendly(self, x): 
47          if self.n == 2: 
48              theta = np.arctan2(x[1, 0], x[0, 0]) 
49              return 'Rot(%.1fdeg)' % np.degrees(theta) 
50          elif self.n == 3: 
51              axis, angle = axis_angle_from_rotation(x) 
52              axisf = S2.friendly(axis) 
53              return 'Rot(%.1fdeg, %s)' % (np.degrees(angle), axisf) 
54          else: 
55              raise ValueError('Not implemented for n>=4.') 
56   
57      @contract(returns='list(belongs)') 
58 -    def interesting_points(self): 
59          points = [] 
60          points.append(self.identity()) 
61          if self.n == 2: 
62              points.append(rot2d(np.pi)) 
63              points.append(rot2d(np.pi / 2)) 
64              points.append(rot2d(-np.pi / 3)) 
65          if self.n == 3: 
66              points.append(rotation_from_axis_angle(np.array([0, 0, 1]), 
67                                                     np.pi / 2)) 
68              points.append(rotation_from_axis_angle(np.array([0, 0, 1]), np.pi)) 
69              points.append(rotation_from_axis_angle(np.array([0, 1, 0]), 
70                                                     np.pi / 2)) 
71              points.append(rotation_from_axis_angle(np.array([0, 1, 0]), np.pi)) 
72              points.append(rotation_from_axis_angle(np.array([1, 0, 0]), 
73                                                     np.pi / 2)) 
74              points.append(rotation_from_axis_angle(np.array([1, 0, 0]), np.pi)) 
75   
76          return points 
77