Source Code for Module geometry.spheres

  1  from . import (safe_arccos, normalize_length, contract, new_contract, 
  2       np) 
  3   
  4  from .utils import assert_allclose 
  5   
  6   
  7  @new_contract 
  8  @contract(x='array[N],N>0') 
  9 -def unit_length(x): 
 10      ''' Checks that the value is a 1D vector with unit length in the 2 norm.''' 
 11      assert_allclose(1, np.linalg.norm(x), rtol=1e-5) # XXX: 
 12   
 13  new_contract('direction', 'array[3], unit_length') 
 14   
 15  new_contract('S1', 'array[2],unit_length') 
 16  new_contract('S2', 'array[3],unit_length') 
 17   
 18   
 19  @new_contract 
 20  @contract(X='array[KxN],K>0,N>0') 
 21 -def directions(X): 
 22      ''' Checks that every column has unit length. ''' 
 23      norm = (X * X).sum(axis=0) 
 24      assert_allclose(1, norm, rtol=1e-5) # XXX: 
 25   
 26   
 27  @contract(s='array[K],K>=2', v='array[K]') 
 28 -def assert_orthogonal(s, v): 
 29      ''' Checks that two vectors are orthogonal. ''' 
 30      dot = (v * s).sum() 
 31      if not np.allclose(dot, 0): 
 32          angle = np.arccos(dot / (np.linalg.norm(v) * np.linalg.norm(s))) 
 33          msg = ('Angle is %.2f deg between %s and %s.' 
 34                 % (np.degrees(angle), s, v)) 
 35          assert_allclose(dot, 0, err_msg=msg) 
 36   
 37   
 38  @contract(x='array[N]', returns='array[N](>=-pi,<pi)') 
 39 -def normalize_pi(x): 
 40      ''' Normalizes the entries in *x* in the interval :math:`[-pi,pi)`. ''' 
 41      angle = np.arctan2(np.sin(x), np.cos(x)) # in [-pi, pi] 
 42      angle[angle == np.pi] = -np.pi 
 43      return angle 
 44   
 45   
 46  @contract(x='float', returns='>=-pi,<pi') 
 47 -def normalize_pi_scalar(x): # TODO: is this the best solution 
 48      angle = np.arctan2(np.sin(x), np.cos(x)) # in [-pi, pi] 
 49      if angle == np.pi: 
 50          return -np.pi 
 51      return angle 
 52   
 53   
 54  @contract(returns='direction') 
 55 -def default_axis(): 
 56      '''  
 57          Returns the axis to use when any will do.  
 58           
 59          For example, the identity is represented by 
 60          a rotation of 0 degrees around *any* axis. If an *(axis,angle)* 
 61          representation is requested, the axis will be given by 
 62          *default_axis()*.  
 63      ''' 
 64      return np.array([0.0, 0.0, 1.0]) 
 65   
 66   
 67  @contract(returns='direction') 
 68 -def default_axis_orthogonal(): 
 69      '''  
 70          Returns an axis orthogonal to the one returned  
 71          by :py:func:`default_axis`.  
 72           
 73          Use this when you need a couple of arbitrary orthogonal axes. 
 74      ''' 
 75      return np.array([0.0, 1.0, 0.0]) 
 76   
 77   
 78  @contract(s1='array[K],unit_length', 
 79             s2='array[K],unit_length', returns='float,>=0,<=pi') 
 80 -def geodesic_distance_on_sphere(s1, s2): 
 81      ''' Returns the geodesic distance between two points on the sphere. ''' 
 82      # special case: return a 0 (no precision issues)  
 83      # if the vectors are the same 
 84      if (s1 == s2).all(): 
 85          return 0.0 
 86      dot_product = (s1 * s2).sum() 
 87      return safe_arccos(dot_product) 
 88   
 89   
 90  @contract(S='directions', returns='float,>=0,<=pi') 
 91 -def distribution_radius(S): 
 92      '''  
 93          Returns the radius of the given directions distribution. 
 94           
 95          The radius is defined as the minimum *r* such that there exists a  
 96          point *s* in *S* such that all distances are within *r* from *s*.  
 97           
 98          .. math:: \\textsf{radius} = \\min \\{ r | \\exists s :  
 99                    \\forall x \\in S : d(s,x) <= r \\} 
100      ''' 
101      D = np.arccos(np.clip(np.dot(S.T, S), -1, 1)) # XXX: repeated 
102      distances = D.max(axis=0) 
103      center = np.argmin(distances) 
104      return distances[center] 
105   
106   
107  @contract(S='array[3xK],directions', s='direction', 
108             returns='array[K](>=0,<=pi)') 
109 -def distances_from(S, s): 
110      '''  
111          Returns the geodesic distances on the sphere from a set of 
112          points *S* to a given point *s*.  
113           
114      ''' 
115      return np.arccos(np.clip(np.dot(s, S), -1, 1)) # XXX: repeated 
116   
117   
118  @contract(ndim='(2|3),K', returns='array[K],unit_length') 
119 -def random_direction(ndim=3): 
120      ''' 
121          Generates a random direction in :math:`\\sphere^{\\ndim-1}`.  
122       
123          Currently only implemented for 2D and 3D. 
124      ''' 
125      if ndim == 3: 
126          z = np.random.uniform(-1, +1) 
127          t = np.random.uniform(0, 2 * np.pi) 
128          r = np.sqrt(1 - z ** 2) 
129          x = r * np.cos(t) 
130          y = r * np.sin(t) 
131          return np.array([x, y, z]) 
132      elif ndim == 2: 
133          theta = np.random.uniform(0, 2 * np.pi) 
134          return np.array([np.cos(theta), np.sin(theta)]) 
135      else: 
136          assert False, 'Not implemented' 
137   
138   
139  @contract(N='int,>0,N', ndim="2|3", returns='array[3xN]') 
140 -def random_directions(N, ndim=3): 
141      ''' Returns a set of random directions. ''' 
142      return np.vstack([random_direction(ndim) for _ in range(N)]).T 
143   
144   
145  @contract(s='direction', returns='direction') 
146 -def any_distant_direction(s): 
147      ''' Returns a direction distant from both *s* and *-s*. ''' 
148      z = default_axis() 
149      d = geodesic_distance_on_sphere(s, z) 
150      # TODO: make this a global parameter 
151      limit = 1.0 / 6.0 * np.pi 
152      if min(d, np.pi - d) < limit: 
153          z = default_axis_orthogonal() 
154      return z 
155   
156   
157  @contract(s='direction', returns='direction') 
158 -def any_orthogonal_direction(s): 
159      ''' Returns any axis orthogonal to *s* (not necessarily random). ''' 
160      # choose a vector far away 
161      z = any_distant_direction(s) 
162      # z ^ s is orthogonal to s 
163      x = np.cross(z, s) 
164      v = x / np.linalg.norm(x) 
165      return v 
166   
167   
168  @contract(s='array[K],unit_length,(K=2|K=3)', returns='array[K],unit_length') 
169 -def random_orthogonal_direction(s): 
170      '''  
171          Returns a random axis orthogonal to *s*  
172          (only implemented for circle and sphere).  
173      ''' 
174      from .rotations import rot2d, rotation_from_axis_angle 
175   
176      if s.size == 2: 
177          theta = np.sign(np.random.uniform() - 0.5) * np.pi / 2 
178          return np.dot(rot2d(theta), s) 
179      elif s.size == 3: 
180          # get any axis orthogonal to s 
181          z = any_orthogonal_direction(s) 
182          # rotate this axis around s by a random amount 
183          angle = np.random.uniform(0, 2 * np.pi) 
184          R = rotation_from_axis_angle(s, angle) 
185          z2 = np.dot(R, z) 
186          return z2 
187      else: 
188          assert False, 'Not implemented' 
189   
190   
191  @contract(s1='array[K],unit_length', s2='array[K],unit_length', 
192            t='number,>=0,<=1') 
193 -def slerp(s1, s2, t): 
194      ''' Spherical interpolation between two points on a hypersphere. ''' 
195      omega = np.arccos(np.dot(s1 / np.linalg.norm(s1), s2 / np.linalg.norm(s2))) 
196      so = np.sin(omega) 
197      if np.abs(so) < 1e-18: # XXX thresholds 
198          return s1 
199      return np.sin((1.0 - t) * omega) / so * s1 + np.sin(t * omega) / so * s2 
200   
201   
202  @contract(ndim='(2|3),K', 
203             radius='number,>0,<=pi', 
204             num_points='int,>0', 
205             center='None|(array[K],unit_length)', 
206             returns='array[KxN],directions') 
207 -def random_directions_bounded(ndim, radius, num_points, center=None): 
208      ''' 
209          Returns a random distribution of points in :math:`\\sphere^{\\ndim-1}`. 
210          within a certain radius from the point *center*.  
211           
212          The points will be distributed uniformly in that area of the sphere. 
213          If *center* is not passed, it will be a random direction.  
214      ''' 
215      from .rotations import rot2d, rotation_from_axis_angle 
216   
217      if center is None: 
218          center = random_direction(ndim) 
219   
220      directions = np.empty((ndim, num_points)) 
221      for i in range(num_points): 
222          # move the center of a random amount  
223          if ndim == 3: 
224              # sample axis orthogonal to the center 
225              axis = random_orthogonal_direction(center) 
226              z = np.random.uniform(0, 1) * spherical_cap_area(radius) 
227              distance = spherical_cap_with_area(z) 
228              R = rotation_from_axis_angle(axis, distance) 
229          elif ndim == 2: 
230              angle = np.random.uniform(-radius, radius) 
231              R = rot2d(angle) 
232          else: 
233              assert False 
234          direction = np.dot(R, center) 
235          directions[:, i] = direction 
236   
237      return sorted_directions(directions) 
238   
239   
240  @contract(S='array[KxN],(K=2|K=3),directions', 
241            returns='array[KxN], directions') 
242 -def sorted_directions(S, num_around=15): 
243      '''  
244          Rearranges the directions in *S* in a better order for visualization. 
245       
246          In 2D, sorts the directions using their angle.  
247           
248          In 3D, it tries to do a pleasant elicoidal arrangement  
249          with **num_around** spires. 
250           
251      ''' 
252      if S.shape[0] == 2: 
253          # XXX check nonzero 
254          center = np.arctan2(S[1, :].sum(), S[0, :].sum()) 
255          angles = np.arctan2(S[1, :], S[0, :]) 
256          diffs = normalize_pi(angles - center) 
257          sorted_d = np.sort(diffs) 
258          final = center + sorted_d 
259          return np.vstack((np.cos(final), np.sin(final))) 
260      else: 
261          # find center of distribution 
262          center = normalize_length(S.sum(axis=1)) 
263          # compute the distances from the center 
264          distance = distances_from(S, center) 
265          # compute the phase from an arbitrary axis 
266          axis = any_orthogonal_direction(center) 
267          phase = distances_from(S, axis) 
268          # normalize distances and phase in [0,1] 
269          phase = (normalize_pi(phase) + np.pi) / (2 * np.pi) 
270          distance /= distance.max() 
271   
272          score = distance * num_around + phase 
273   
274          order = np.argsort(score) 
275   
276          ordered = S[:, order] 
277          return ordered 
278   
279   
280 -def sphere_area(r=1): 
281      ''' Returns the area of a sphere of the given radius. ''' 
282      return 4 * np.pi * (r ** 2) 
283   
284   
285 -def spherical_cap_area(cap_radius): 
286      '''  
287          Returns the area of a spherical cap on the unit sphere  
288          of the given radius.  
289       
290          See figure at http://mathworld.wolfram.com/SphericalCap.html 
291      ''' 
292      h = 1 - np.cos(cap_radius) 
293      a = np.sin(cap_radius) 
294      A = np.pi * (a ** 2 + h ** 2) 
295      return A 
296   
297   
298 -def spherical_cap_with_area(cap_area): 
299      '''  
300          Returns the radius of a spherical cap of the given area.  
301       
302          See http://www.springerlink.com/content/3521h167300g7v62/ 
303      ''' 
304      A = cap_area 
305      L = np.sqrt(A / np.pi) 
306      h = L ** 2 / 2 
307      r = np.arccos(1 - h) 
308      return r 
309   
310   
311  @contract(S='array[KxN],K>=2', returns='array[KxN]') 
312 -def project_vectors_onto_sphere(S, atol=1e-7): 
313      K, N = S.shape 
314      coords_proj = np.zeros((K, N)) 
315      for i in range(N): 
316          v = S[:, i] 
317          nv = np.linalg.norm(v) 
318          if np.fabs(nv) < atol: 
319              raise ValueError('Vector too small: %s' % v) 
320          coords_proj[:, i] = v / nv 
321      return coords_proj 
322