Source Code for Module geometry.manifolds.sphere

  1  from . import DifferentiableManifold, np, contract, check 
  2  from .. import (any_orthogonal_direction, geodesic_distance_on_sphere, 
  3      random_direction, normalize_length, normalize_length_or_zero, 
  4      rotation_from_axis_angle, rot2d) 
  5  from numpy.core.numeric import outer 
  6   
  7   
  8 -class Sphere(DifferentiableManifold): 
  9      ''' These are hyperspheres of unit radius. ''' 
 10   
 11      norm_rtol = 1e-5 
 12      atol_geodesic_distance = 1e-8 
 13   
 14      @contract(order='(1|2)') 
 15 -    def __init__(self, order): 
 16          DifferentiableManifold.__init__(self, dimension=order) 
 17          self.N = order + 1 
 18   
 19 -    def __repr__(self): 
 20          return 'S%s' % (self.dimension) 
 21   
 22      @contract(a='belongs', b='belongs', returns='>=0') 
 23 -    def distance(self, a, b): 
 24          return geodesic_distance_on_sphere(a, b) 
 25   
 26      @contract(base='belongs', target='belongs', returns='belongs_ts') 
 27 -    def logmap(self, base, target): 
 28          # TODO: create S1_logmap(base, target) 
 29          # XXX: what should we do in the case there is more than one logmap? 
 30          d = geodesic_distance_on_sphere(target, base) 
 31          if np.allclose(d, np.pi, atol=self.atol_geodesic_distance): 
 32              if self.N == 2: 
 33                  xp = np.array([base[1], -base[0]]) 
 34              elif self.N == 3: 
 35                  xp = any_orthogonal_direction(base) 
 36          else: 
 37              x = target - base 
 38              base, xp = self.project_ts((base, x)) 
 39          xp = normalize_length_or_zero(xp) 
 40          xp *= geodesic_distance_on_sphere(target, base) 
 41          return base, xp 
 42   
 43      @contract(bv='belongs_ts', returns='belongs') 
 44 -    def expmap(self, bv): 
 45          base, vel = bv 
 46          angle = np.linalg.norm(vel) 
 47          if angle == 0: # XXX: tolerance 
 48              return base 
 49   
 50          if self.N == 2: 
 51              direction = -np.sign(vel[0] * base[1] - base[0] * vel[1]) 
 52              R = rot2d(angle * direction) 
 53              result = np.dot(R, base) 
 54          elif self.N == 3: 
 55              axis = -np.cross(normalize_length(vel), base) 
 56              axis = normalize_length(axis) 
 57              R = rotation_from_axis_angle(axis, angle) 
 58              result = np.dot(R, base) 
 59          else: 
 60              assert False 
 61   
 62          return result 
 63   
 64      @contract(bve='tuple(belongs, *)') 
 65 -    def project_ts(self, bve): # TODO: test 
 66          base, vel = bve 
 67          P = np.eye(self.N) - outer(base, base) 
 68          return base, np.dot(P, vel) 
 69   
 70 -    def belongs(self, x): 
 71          check('array[N],unit_length', x) 
 72   
 73      @contract(returns='belongs') 
 74 -    def sample_uniform(self): 
 75          return random_direction(self.N) 
 76   
 77      @contract(returns='list(belongs)') 
 78 -    def interesting_points(self): 
 79          if self.N == 2: 
 80              points = [] 
 81              points.append(np.array([-1, 0])) 
 82              points.append(np.array([0, 1])) 
 83              points.append(np.array([0, -1])) 
 84              points.append(np.array([+1, 0])) 
 85              points.append(np.array([np.sqrt(2) / 2, np.sqrt(2) / 2])) 
 86              return points 
 87          elif self.N == 3: 
 88              points = [] 
 89              points.append(np.array([-1, 0, 0])) 
 90              points.append(np.array([0, 1, 0])) 
 91              points.append(np.array([0, 0, 1])) 
 92              points.append(np.array([0, 0, -1])) 
 93              points.append(np.array([np.sqrt(2) / 2, np.sqrt(2) / 2, 0])) 
 94              return points 
 95          else: 
 96              assert False 
 97   
 98 -    def friendly(self, x): 
 99          if self.N == 2: 
100              theta = np.arctan2(x[1], x[0]) 
101              return 'Dir(%6.1fdeg)' % np.degrees(theta) 
102          elif self.N == 3: 
103              theta = np.arctan2(x[1], x[0]) 
104              elevation = np.arcsin(x[2]) 
105              return 'Dir(%6.1fdeg,el:%5.1fdeg)' % (np.degrees(theta), 
106                                                        np.degrees(elevation)) 
107          else: 
108              assert False 
109   
110   
111 -class Sphere1(DifferentiableManifold): 
112   
113      norm_rtol = 1e-5 
114      atol_geodesic_distance = 1e-8 
115   
116 -    def __init__(self): 
117          DifferentiableManifold.__init__(self, dimension=1) 
118   
119 -    def __repr__(self): 
120          return 'S1' 
121   
122      @contract(a='S1', b='S1', returns='>=0') 
123 -    def distance(self, a, b): 
124          return geodesic_distance_on_sphere(a, b) 
125   
126      @contract(base='S1', target='S1', returns='belongs_ts') 
127 -    def logmap(self, base, target): 
128          # TODO: create S1_logmap(base, target) 
129          # XXX: what should we do in the case there is more than one logmap? 
130          d = geodesic_distance_on_sphere(target, base) 
131          if np.allclose(d, np.pi, atol=self.atol_geodesic_distance): 
132              xp = np.array([base[1], -base[0]]) 
133          else: 
134              x = target - base 
135              base, xp = self.project_ts((base, x)) 
136          xp = normalize_length_or_zero(xp) 
137          xp *= geodesic_distance_on_sphere(target, base) 
138          return base, xp 
139   
140      @contract(bv='belongs_ts', returns='belongs') 
141 -    def expmap(self, bv): 
142          base, vel = bv 
143          angle = np.linalg.norm(vel) 
144          if angle == 0: # XXX: tolerance 
145              return base 
146          direction = -np.sign(vel[0] * base[1] - base[0] * vel[1]) 
147          R = rot2d(angle * direction) 
148          result = np.dot(R, base) 
149          return result 
150   
151      @contract(bx='tuple(belongs, *)') 
152 -    def project_ts(self, bx): # TODO: test 
153          base, x = bx 
154          P = np.eye(2) - outer(base, base) 
155          return base, np.dot(P, x) 
156   
157      @contract(x='S1') 
158 -    def belongs(self, x): 
159          pass 
160   
161 -    def sample_uniform(self): 
162          return random_direction(2) 
163   
164 -    def interesting_points(self): 
165          points = [] 
166          points.append(np.array([-1, 0])) 
167          points.append(np.array([0, 1])) 
168          points.append(np.array([0, -1])) 
169          points.append(np.array([+1, 0])) 
170          points.append(np.array([np.sqrt(2) / 2, np.sqrt(2) / 2])) 
171          return points 
172   
173 -    def friendly(self, x): 
174          theta = np.arctan2(x[1], x[0]) 
175          return 'Dir(%6.1fdeg)' % np.degrees(theta) 
176