Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/spatial/_geometric_slerp.py : 15%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1from __future__ import division, print_function, absolute_import
3__all__ = ['geometric_slerp']
5import warnings
7import numpy as np
8from scipy.spatial.distance import euclidean
11def _geometric_slerp(start, end, t):
12 # create an orthogonal basis using QR decomposition
13 basis = np.vstack([start, end])
14 Q, R = np.linalg.qr(basis.T)
15 signs = 2 * (np.diag(R) >= 0) - 1
16 Q = Q.T * signs.T[:, np.newaxis]
17 R = R.T * signs.T[:, np.newaxis]
19 # calculate the angle between `start` and `end`
20 c = np.dot(start, end)
21 s = np.linalg.det(R)
22 omega = np.arctan2(s, c)
24 # interpolate
25 start, end = Q
26 s = np.sin(t * omega)
27 c = np.cos(t * omega)
28 return start * c[:, np.newaxis] + end * s[:, np.newaxis]
31def geometric_slerp(start,
32 end,
33 t,
34 tol=1e-7):
35 """
36 Geometric spherical linear interpolation.
38 The interpolation occurs along a unit-radius
39 great circle arc in arbitrary dimensional space.
41 Parameters
42 ----------
43 start : (n_dimensions, ) array-like
44 Single n-dimensional input coordinate in a 1-D array-like
45 object. `n` must be greater than 1.
46 end : (n_dimensions, ) array-like
47 Single n-dimensional input coordinate in a 1-D array-like
48 object. `n` must be greater than 1.
49 t: float or (n_points,) array-like
50 A float or array-like of doubles representing interpolation
51 parameters, with values required in the inclusive interval
52 between 0 and 1. A common approach is to generate the array
53 with ``np.linspace(0, 1, n_pts)`` for linearly spaced points.
54 Ascending, descending, and scrambled orders are permitted.
55 tol: float
56 The absolute tolerance for determining if the start and end
57 coordinates are antipodes.
59 Returns
60 -------
61 result : (t.size, D)
62 An array of doubles containing the interpolated
63 spherical path and including start and
64 end when 0 and 1 t are used. The
65 interpolated values should correspond to the
66 same sort order provided in the t array. The result
67 may be 1-dimensional if ``t`` is a float.
69 Raises
70 ------
71 ValueError
72 If ``start`` and ``end`` are antipodes, not on the
73 unit n-sphere, or for a variety of degenerate conditions.
75 Notes
76 -----
77 The implementation is based on the mathematical formula provided in [1]_,
78 and the first known presentation of this algorithm, derived from study of
79 4-D geometry, is credited to Glenn Davis in a footnote of the original
80 quaternion Slerp publication by Ken Shoemake [2]_.
82 .. versionadded:: 1.5.0
84 References
85 ----------
86 .. [1] https://en.wikipedia.org/wiki/Slerp#Geometric_Slerp
87 .. [2] Ken Shoemake (1985) Animating rotation with quaternion curves.
88 ACM SIGGRAPH Computer Graphics, 19(3): 245-254.
90 See Also
91 --------
92 scipy.spatial.transform.Slerp : 3-D Slerp that works with quaternions
94 Examples
95 --------
96 Interpolate four linearly-spaced values on the circumference of
97 a circle spanning 90 degrees:
99 >>> from scipy.spatial import geometric_slerp
100 >>> import matplotlib.pyplot as plt
101 >>> fig = plt.figure()
102 >>> ax = fig.add_subplot(111)
103 >>> start = np.array([1, 0])
104 >>> end = np.array([0, 1])
105 >>> t_vals = np.linspace(0, 1, 4)
106 >>> result = geometric_slerp(start,
107 ... end,
108 ... t_vals)
110 The interpolated results should be at 30 degree intervals
111 recognizable on the unit circle:
113 >>> ax.scatter(result[...,0], result[...,1], c='k')
114 >>> circle = plt.Circle((0, 0), 1, color='grey')
115 >>> ax.add_artist(circle)
116 >>> ax.set_aspect('equal')
117 >>> plt.show()
119 Attempting to interpolate between antipodes on a circle is
120 ambiguous because there are two possible paths, and on a
121 sphere there are infinite possible paths on the geodesic surface.
122 Nonetheless, one of the ambiguous paths is returned along
123 with a warning:
125 >>> opposite_pole = np.array([-1, 0])
126 >>> with np.testing.suppress_warnings() as sup:
127 ... sup.filter(UserWarning)
128 ... geometric_slerp(start,
129 ... opposite_pole,
130 ... t_vals)
131 array([[ 1.00000000e+00, 0.00000000e+00],
132 [ 5.00000000e-01, 8.66025404e-01],
133 [-5.00000000e-01, 8.66025404e-01],
134 [-1.00000000e+00, 1.22464680e-16]])
136 Extend the original example to a sphere and plot interpolation
137 points in 3D:
139 >>> from mpl_toolkits.mplot3d import proj3d
140 >>> fig = plt.figure()
141 >>> ax = fig.add_subplot(111, projection='3d')
143 Plot the unit sphere for reference (optional):
145 >>> u = np.linspace(0, 2 * np.pi, 100)
146 >>> v = np.linspace(0, np.pi, 100)
147 >>> x = np.outer(np.cos(u), np.sin(v))
148 >>> y = np.outer(np.sin(u), np.sin(v))
149 >>> z = np.outer(np.ones(np.size(u)), np.cos(v))
150 >>> ax.plot_surface(x, y, z, color='y', alpha=0.1)
152 Interpolating over a larger number of points
153 may provide the appearance of a smooth curve on
154 the surface of the sphere, which is also useful
155 for discretized integration calculations on a
156 sphere surface:
158 >>> start = np.array([1, 0, 0])
159 >>> end = np.array([0, 0, 1])
160 >>> t_vals = np.linspace(0, 1, 200)
161 >>> result = geometric_slerp(start,
162 ... end,
163 ... t_vals)
164 >>> ax.plot(result[...,0],
165 ... result[...,1],
166 ... result[...,2],
167 ... c='k')
168 >>> plt.show()
169 """
171 start = np.asarray(start, dtype=np.float64)
172 end = np.asarray(end, dtype=np.float64)
174 if start.ndim != 1 or end.ndim != 1:
175 raise ValueError("Start and end coordinates "
176 "must be one-dimensional")
178 if start.size != end.size:
179 raise ValueError("The dimensions of start and "
180 "end must match (have same size)")
182 if start.size < 2 or end.size < 2:
183 raise ValueError("The start and end coordinates must "
184 "both be in at least two-dimensional "
185 "space")
187 if np.array_equal(start, end):
188 return [start] * np.asarray(t).size
190 # for points that violate equation for n-sphere
191 for coord in [start, end]:
192 if not np.allclose(np.linalg.norm(coord), 1.0,
193 rtol=1e-9,
194 atol=0):
195 raise ValueError("start and end are not"
196 " on a unit n-sphere")
198 if not isinstance(tol, float):
199 raise ValueError("tol must be a float")
200 else:
201 tol = np.fabs(tol)
203 coord_dist = euclidean(start, end)
205 # diameter of 2 within tolerance means antipodes, which is a problem
206 # for all unit n-spheres (even the 0-sphere would have an ambiguous path)
207 if np.allclose(coord_dist, 2.0, rtol=0, atol=tol):
208 warnings.warn("start and end are antipodes"
209 " using the specified tolerance;"
210 " this may cause ambiguous slerp paths")
212 t = np.asarray(t, dtype=np.float64)
214 if t.size == 0:
215 return np.empty((0, start.size))
217 if t.min() < 0 or t.max() > 1:
218 raise ValueError("interpolation parameter must be in [0, 1]")
220 if t.ndim == 0:
221 return _geometric_slerp(start,
222 end,
223 np.atleast_1d(t)).ravel()
224 else:
225 return _geometric_slerp(start,
226 end,
227 t)