Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/spatial/_geometric_slerp.py: 17%

1from __future__ import annotations 

2 

3__all__ = ['geometric_slerp'] 

4 

5import warnings 

6from typing import TYPE_CHECKING 

7 

8import numpy as np 

9from scipy.spatial.distance import euclidean 

10 

11if TYPE_CHECKING: 

12 import numpy.typing as npt 

13 

14 

15def _geometric_slerp(start, end, t): 

16 # create an orthogonal basis using QR decomposition 

17 basis = np.vstack([start, end]) 

18 Q, R = np.linalg.qr(basis.T) 

19 signs = 2 * (np.diag(R) >= 0) - 1 

20 Q = Q.T * signs.T[:, np.newaxis] 

21 R = R.T * signs.T[:, np.newaxis] 

22 

23 # calculate the angle between `start` and `end` 

24 c = np.dot(start, end) 

25 s = np.linalg.det(R) 

26 omega = np.arctan2(s, c) 

27 

28 # interpolate 

29 start, end = Q 

30 s = np.sin(t * omega) 

31 c = np.cos(t * omega) 

32 return start * c[:, np.newaxis] + end * s[:, np.newaxis] 

33 

34 

35def geometric_slerp( 

36 start: npt.ArrayLike, 

37 end: npt.ArrayLike, 

38 t: npt.ArrayLike, 

39 tol: float = 1e-7, 

40) -> np.ndarray: 

41 """ 

42 Geometric spherical linear interpolation. 

43 

44 The interpolation occurs along a unit-radius 

45 great circle arc in arbitrary dimensional space. 

46 

47 Parameters 

48 ---------- 

49 start : (n_dimensions, ) array-like 

50 Single n-dimensional input coordinate in a 1-D array-like 

51 object. `n` must be greater than 1. 

52 end : (n_dimensions, ) array-like 

53 Single n-dimensional input coordinate in a 1-D array-like 

54 object. `n` must be greater than 1. 

55 t : float or (n_points,) 1D array-like 

56 A float or 1D array-like of doubles representing interpolation 

57 parameters, with values required in the inclusive interval 

58 between 0 and 1. A common approach is to generate the array 

59 with ``np.linspace(0, 1, n_pts)`` for linearly spaced points. 

60 Ascending, descending, and scrambled orders are permitted. 

61 tol : float 

62 The absolute tolerance for determining if the start and end 

63 coordinates are antipodes. 

64 

65 Returns 

66 ------- 

67 result : (t.size, D) 

68 An array of doubles containing the interpolated 

69 spherical path and including start and 

70 end when 0 and 1 t are used. The 

71 interpolated values should correspond to the 

72 same sort order provided in the t array. The result 

73 may be 1-dimensional if ``t`` is a float. 

74 

75 Raises 

76 ------ 

77 ValueError 

78 If ``start`` and ``end`` are antipodes, not on the 

79 unit n-sphere, or for a variety of degenerate conditions. 

80 

81 See Also 

82 -------- 

83 scipy.spatial.transform.Slerp : 3-D Slerp that works with quaternions 

84 

85 Notes 

86 ----- 

87 The implementation is based on the mathematical formula provided in [1]_, 

88 and the first known presentation of this algorithm, derived from study of 

89 4-D geometry, is credited to Glenn Davis in a footnote of the original 

90 quaternion Slerp publication by Ken Shoemake [2]_. 

91 

92 .. versionadded:: 1.5.0 

93 

94 References 

95 ---------- 

96 .. [1] https://en.wikipedia.org/wiki/Slerp#Geometric_Slerp 

97 .. [2] Ken Shoemake (1985) Animating rotation with quaternion curves. 

98 ACM SIGGRAPH Computer Graphics, 19(3): 245-254. 

99 

100 Examples 

101 -------- 

102 Interpolate four linearly-spaced values on the circumference of 

103 a circle spanning 90 degrees: 

104 

105 >>> import numpy as np 

106 >>> from scipy.spatial import geometric_slerp 

107 >>> import matplotlib.pyplot as plt 

108 >>> fig = plt.figure() 

109 >>> ax = fig.add_subplot(111) 

110 >>> start = np.array([1, 0]) 

111 >>> end = np.array([0, 1]) 

112 >>> t_vals = np.linspace(0, 1, 4) 

113 >>> result = geometric_slerp(start, 

114 ... end, 

115 ... t_vals) 

116 

117 The interpolated results should be at 30 degree intervals 

118 recognizable on the unit circle: 

119 

120 >>> ax.scatter(result[...,0], result[...,1], c='k') 

121 >>> circle = plt.Circle((0, 0), 1, color='grey') 

122 >>> ax.add_artist(circle) 

123 >>> ax.set_aspect('equal') 

124 >>> plt.show() 

125 

126 Attempting to interpolate between antipodes on a circle is 

127 ambiguous because there are two possible paths, and on a 

128 sphere there are infinite possible paths on the geodesic surface. 

129 Nonetheless, one of the ambiguous paths is returned along 

130 with a warning: 

131 

132 >>> opposite_pole = np.array([-1, 0]) 

133 >>> with np.testing.suppress_warnings() as sup: 

134 ... sup.filter(UserWarning) 

135 ... geometric_slerp(start, 

136 ... opposite_pole, 

137 ... t_vals) 

138 array([[ 1.00000000e+00, 0.00000000e+00], 

139 [ 5.00000000e-01, 8.66025404e-01], 

140 [-5.00000000e-01, 8.66025404e-01], 

141 [-1.00000000e+00, 1.22464680e-16]]) 

142 

143 Extend the original example to a sphere and plot interpolation 

144 points in 3D: 

145 

146 >>> from mpl_toolkits.mplot3d import proj3d 

147 >>> fig = plt.figure() 

148 >>> ax = fig.add_subplot(111, projection='3d') 

149 

150 Plot the unit sphere for reference (optional): 

151 

152 >>> u = np.linspace(0, 2 * np.pi, 100) 

153 >>> v = np.linspace(0, np.pi, 100) 

154 >>> x = np.outer(np.cos(u), np.sin(v)) 

155 >>> y = np.outer(np.sin(u), np.sin(v)) 

156 >>> z = np.outer(np.ones(np.size(u)), np.cos(v)) 

157 >>> ax.plot_surface(x, y, z, color='y', alpha=0.1) 

158 

159 Interpolating over a larger number of points 

160 may provide the appearance of a smooth curve on 

161 the surface of the sphere, which is also useful 

162 for discretized integration calculations on a 

163 sphere surface: 

164 

165 >>> start = np.array([1, 0, 0]) 

166 >>> end = np.array([0, 0, 1]) 

167 >>> t_vals = np.linspace(0, 1, 200) 

168 >>> result = geometric_slerp(start, 

169 ... end, 

170 ... t_vals) 

171 >>> ax.plot(result[...,0], 

172 ... result[...,1], 

173 ... result[...,2], 

174 ... c='k') 

175 >>> plt.show() 

176 """ 

177 

178 start = np.asarray(start, dtype=np.float64) 

179 end = np.asarray(end, dtype=np.float64) 

180 t = np.asarray(t) 

181 

182 if t.ndim > 1: 

183 raise ValueError("The interpolation parameter " 

184 "value must be one dimensional.") 

185 

186 if start.ndim != 1 or end.ndim != 1: 

187 raise ValueError("Start and end coordinates " 

188 "must be one-dimensional") 

189 

190 if start.size != end.size: 

191 raise ValueError("The dimensions of start and " 

192 "end must match (have same size)") 

193 

194 if start.size < 2 or end.size < 2: 

195 raise ValueError("The start and end coordinates must " 

196 "both be in at least two-dimensional " 

197 "space") 

198 

199 if np.array_equal(start, end): 

200 return np.linspace(start, start, t.size) 

201 

202 # for points that violate equation for n-sphere 

203 for coord in [start, end]: 

204 if not np.allclose(np.linalg.norm(coord), 1.0, 

205 rtol=1e-9, 

206 atol=0): 

207 raise ValueError("start and end are not" 

208 " on a unit n-sphere") 

209 

210 if not isinstance(tol, float): 

211 raise ValueError("tol must be a float") 

212 else: 

213 tol = np.fabs(tol) 

214 

215 coord_dist = euclidean(start, end) 

216 

217 # diameter of 2 within tolerance means antipodes, which is a problem 

218 # for all unit n-spheres (even the 0-sphere would have an ambiguous path) 

219 if np.allclose(coord_dist, 2.0, rtol=0, atol=tol): 

220 warnings.warn("start and end are antipodes" 

221 " using the specified tolerance;" 

222 " this may cause ambiguous slerp paths") 

223 

224 t = np.asarray(t, dtype=np.float64) 

225 

226 if t.size == 0: 

227 return np.empty((0, start.size)) 

228 

229 if t.min() < 0 or t.max() > 1: 

230 raise ValueError("interpolation parameter must be in [0, 1]") 

231 

232 if t.ndim == 0: 

233 return _geometric_slerp(start, 

234 end, 

235 np.atleast_1d(t)).ravel() 

236 else: 

237 return _geometric_slerp(start, 

238 end, 

239 t)