proj3d.py 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259
  1. """
  2. Various transforms used for by the 3D code
  3. """
  4. import numpy as np
  5. from matplotlib import _api
  6. def world_transformation(xmin, xmax,
  7. ymin, ymax,
  8. zmin, zmax, pb_aspect=None):
  9. """
  10. Produce a matrix that scales homogeneous coords in the specified ranges
  11. to [0, 1], or [0, pb_aspect[i]] if the plotbox aspect ratio is specified.
  12. """
  13. dx = xmax - xmin
  14. dy = ymax - ymin
  15. dz = zmax - zmin
  16. if pb_aspect is not None:
  17. ax, ay, az = pb_aspect
  18. dx /= ax
  19. dy /= ay
  20. dz /= az
  21. return np.array([[1/dx, 0, 0, -xmin/dx],
  22. [0, 1/dy, 0, -ymin/dy],
  23. [0, 0, 1/dz, -zmin/dz],
  24. [0, 0, 0, 1]])
  25. @_api.deprecated("3.8")
  26. def rotation_about_vector(v, angle):
  27. """
  28. Produce a rotation matrix for an angle in radians about a vector.
  29. """
  30. return _rotation_about_vector(v, angle)
  31. def _rotation_about_vector(v, angle):
  32. """
  33. Produce a rotation matrix for an angle in radians about a vector.
  34. """
  35. vx, vy, vz = v / np.linalg.norm(v)
  36. s = np.sin(angle)
  37. c = np.cos(angle)
  38. t = 2*np.sin(angle/2)**2 # more numerically stable than t = 1-c
  39. R = np.array([
  40. [t*vx*vx + c, t*vx*vy - vz*s, t*vx*vz + vy*s],
  41. [t*vy*vx + vz*s, t*vy*vy + c, t*vy*vz - vx*s],
  42. [t*vz*vx - vy*s, t*vz*vy + vx*s, t*vz*vz + c]])
  43. return R
  44. def _view_axes(E, R, V, roll):
  45. """
  46. Get the unit viewing axes in data coordinates.
  47. Parameters
  48. ----------
  49. E : 3-element numpy array
  50. The coordinates of the eye/camera.
  51. R : 3-element numpy array
  52. The coordinates of the center of the view box.
  53. V : 3-element numpy array
  54. Unit vector in the direction of the vertical axis.
  55. roll : float
  56. The roll angle in radians.
  57. Returns
  58. -------
  59. u : 3-element numpy array
  60. Unit vector pointing towards the right of the screen.
  61. v : 3-element numpy array
  62. Unit vector pointing towards the top of the screen.
  63. w : 3-element numpy array
  64. Unit vector pointing out of the screen.
  65. """
  66. w = (E - R)
  67. w = w/np.linalg.norm(w)
  68. u = np.cross(V, w)
  69. u = u/np.linalg.norm(u)
  70. v = np.cross(w, u) # Will be a unit vector
  71. # Save some computation for the default roll=0
  72. if roll != 0:
  73. # A positive rotation of the camera is a negative rotation of the world
  74. Rroll = _rotation_about_vector(w, -roll)
  75. u = np.dot(Rroll, u)
  76. v = np.dot(Rroll, v)
  77. return u, v, w
  78. def _view_transformation_uvw(u, v, w, E):
  79. """
  80. Return the view transformation matrix.
  81. Parameters
  82. ----------
  83. u : 3-element numpy array
  84. Unit vector pointing towards the right of the screen.
  85. v : 3-element numpy array
  86. Unit vector pointing towards the top of the screen.
  87. w : 3-element numpy array
  88. Unit vector pointing out of the screen.
  89. E : 3-element numpy array
  90. The coordinates of the eye/camera.
  91. """
  92. Mr = np.eye(4)
  93. Mt = np.eye(4)
  94. Mr[:3, :3] = [u, v, w]
  95. Mt[:3, -1] = -E
  96. M = np.dot(Mr, Mt)
  97. return M
  98. @_api.deprecated("3.8")
  99. def view_transformation(E, R, V, roll):
  100. """
  101. Return the view transformation matrix.
  102. Parameters
  103. ----------
  104. E : 3-element numpy array
  105. The coordinates of the eye/camera.
  106. R : 3-element numpy array
  107. The coordinates of the center of the view box.
  108. V : 3-element numpy array
  109. Unit vector in the direction of the vertical axis.
  110. roll : float
  111. The roll angle in radians.
  112. """
  113. u, v, w = _view_axes(E, R, V, roll)
  114. M = _view_transformation_uvw(u, v, w, E)
  115. return M
  116. @_api.deprecated("3.8")
  117. def persp_transformation(zfront, zback, focal_length):
  118. return _persp_transformation(zfront, zback, focal_length)
  119. def _persp_transformation(zfront, zback, focal_length):
  120. e = focal_length
  121. a = 1 # aspect ratio
  122. b = (zfront+zback)/(zfront-zback)
  123. c = -2*(zfront*zback)/(zfront-zback)
  124. proj_matrix = np.array([[e, 0, 0, 0],
  125. [0, e/a, 0, 0],
  126. [0, 0, b, c],
  127. [0, 0, -1, 0]])
  128. return proj_matrix
  129. @_api.deprecated("3.8")
  130. def ortho_transformation(zfront, zback):
  131. return _ortho_transformation(zfront, zback)
  132. def _ortho_transformation(zfront, zback):
  133. # note: w component in the resulting vector will be (zback-zfront), not 1
  134. a = -(zfront + zback)
  135. b = -(zfront - zback)
  136. proj_matrix = np.array([[2, 0, 0, 0],
  137. [0, 2, 0, 0],
  138. [0, 0, -2, 0],
  139. [0, 0, a, b]])
  140. return proj_matrix
  141. def _proj_transform_vec(vec, M):
  142. vecw = np.dot(M, vec)
  143. w = vecw[3]
  144. # clip here..
  145. txs, tys, tzs = vecw[0]/w, vecw[1]/w, vecw[2]/w
  146. return txs, tys, tzs
  147. def _proj_transform_vec_clip(vec, M):
  148. vecw = np.dot(M, vec)
  149. w = vecw[3]
  150. # clip here.
  151. txs, tys, tzs = vecw[0] / w, vecw[1] / w, vecw[2] / w
  152. tis = (0 <= vecw[0]) & (vecw[0] <= 1) & (0 <= vecw[1]) & (vecw[1] <= 1)
  153. if np.any(tis):
  154. tis = vecw[1] < 1
  155. return txs, tys, tzs, tis
  156. def inv_transform(xs, ys, zs, invM):
  157. """
  158. Transform the points by the inverse of the projection matrix, *invM*.
  159. """
  160. vec = _vec_pad_ones(xs, ys, zs)
  161. vecr = np.dot(invM, vec)
  162. if vecr.shape == (4,):
  163. vecr = vecr.reshape((4, 1))
  164. for i in range(vecr.shape[1]):
  165. if vecr[3][i] != 0:
  166. vecr[:, i] = vecr[:, i] / vecr[3][i]
  167. return vecr[0], vecr[1], vecr[2]
  168. def _vec_pad_ones(xs, ys, zs):
  169. return np.array([xs, ys, zs, np.ones_like(xs)])
  170. def proj_transform(xs, ys, zs, M):
  171. """
  172. Transform the points by the projection matrix *M*.
  173. """
  174. vec = _vec_pad_ones(xs, ys, zs)
  175. return _proj_transform_vec(vec, M)
  176. transform = _api.deprecated(
  177. "3.8", obj_type="function", name="transform",
  178. alternative="proj_transform")(proj_transform)
  179. def proj_transform_clip(xs, ys, zs, M):
  180. """
  181. Transform the points by the projection matrix
  182. and return the clipping result
  183. returns txs, tys, tzs, tis
  184. """
  185. vec = _vec_pad_ones(xs, ys, zs)
  186. return _proj_transform_vec_clip(vec, M)
  187. @_api.deprecated("3.8")
  188. def proj_points(points, M):
  189. return _proj_points(points, M)
  190. def _proj_points(points, M):
  191. return np.column_stack(_proj_trans_points(points, M))
  192. @_api.deprecated("3.8")
  193. def proj_trans_points(points, M):
  194. return _proj_trans_points(points, M)
  195. def _proj_trans_points(points, M):
  196. xs, ys, zs = zip(*points)
  197. return proj_transform(xs, ys, zs, M)
  198. @_api.deprecated("3.8")
  199. def rot_x(V, alpha):
  200. cosa, sina = np.cos(alpha), np.sin(alpha)
  201. M1 = np.array([[1, 0, 0, 0],
  202. [0, cosa, -sina, 0],
  203. [0, sina, cosa, 0],
  204. [0, 0, 0, 1]])
  205. return np.dot(M1, V)