proj3d.py 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239
  1. """
  2. Various transforms used for by the 3D code
  3. """
  4. import numpy as np
  5. import numpy.linalg as linalg
  6. from matplotlib import cbook
  7. @cbook.deprecated("3.1")
  8. def line2d(p0, p1):
  9. """
  10. Return 2D equation of line in the form ax+by+c = 0
  11. """
  12. # x + x1 = 0
  13. x0, y0 = p0[:2]
  14. x1, y1 = p1[:2]
  15. #
  16. if x0 == x1:
  17. a = -1
  18. b = 0
  19. c = x1
  20. elif y0 == y1:
  21. a = 0
  22. b = 1
  23. c = -y1
  24. else:
  25. a = y0 - y1
  26. b = x0 - x1
  27. c = x0*y1 - x1*y0
  28. return a, b, c
  29. @cbook.deprecated("3.1")
  30. def line2d_dist(l, p):
  31. """
  32. Distance from line to point
  33. line is a tuple of coefficients a, b, c
  34. """
  35. a, b, c = l
  36. x0, y0 = p
  37. return abs((a*x0 + b*y0 + c) / np.hypot(a, b))
  38. def _line2d_seg_dist(p1, p2, p0):
  39. """distance(s) from line defined by p1 - p2 to point(s) p0
  40. p0[0] = x(s)
  41. p0[1] = y(s)
  42. intersection point p = p1 + u*(p2-p1)
  43. and intersection point lies within segment if u is between 0 and 1
  44. """
  45. x21 = p2[0] - p1[0]
  46. y21 = p2[1] - p1[1]
  47. x01 = np.asarray(p0[0]) - p1[0]
  48. y01 = np.asarray(p0[1]) - p1[1]
  49. u = (x01*x21 + y01*y21) / (x21**2 + y21**2)
  50. u = np.clip(u, 0, 1)
  51. d = np.hypot(x01 - u*x21, y01 - u*y21)
  52. return d
  53. @cbook.deprecated("3.1")
  54. def line2d_seg_dist(p1, p2, p0):
  55. """distance(s) from line defined by p1 - p2 to point(s) p0
  56. p0[0] = x(s)
  57. p0[1] = y(s)
  58. intersection point p = p1 + u*(p2-p1)
  59. and intersection point lies within segment if u is between 0 and 1
  60. """
  61. return _line2d_seg_dist(p1, p2, p0)
  62. @cbook.deprecated("3.1", alternative="np.linalg.norm")
  63. def mod(v):
  64. """3d vector length"""
  65. return np.sqrt(v[0]**2+v[1]**2+v[2]**2)
  66. def world_transformation(xmin, xmax,
  67. ymin, ymax,
  68. zmin, zmax):
  69. dx, dy, dz = (xmax-xmin), (ymax-ymin), (zmax-zmin)
  70. return np.array([[1/dx, 0, 0, -xmin/dx],
  71. [0, 1/dy, 0, -ymin/dy],
  72. [0, 0, 1/dz, -zmin/dz],
  73. [0, 0, 0, 1]])
  74. def view_transformation(E, R, V):
  75. n = (E - R)
  76. ## new
  77. # n /= np.linalg.norm(n)
  78. # u = np.cross(V, n)
  79. # u /= np.linalg.norm(u)
  80. # v = np.cross(n, u)
  81. # Mr = np.diag([1.] * 4)
  82. # Mt = np.diag([1.] * 4)
  83. # Mr[:3,:3] = u, v, n
  84. # Mt[:3,-1] = -E
  85. ## end new
  86. ## old
  87. n = n / np.linalg.norm(n)
  88. u = np.cross(V, n)
  89. u = u / np.linalg.norm(u)
  90. v = np.cross(n, u)
  91. Mr = [[u[0], u[1], u[2], 0],
  92. [v[0], v[1], v[2], 0],
  93. [n[0], n[1], n[2], 0],
  94. [0, 0, 0, 1]]
  95. #
  96. Mt = [[1, 0, 0, -E[0]],
  97. [0, 1, 0, -E[1]],
  98. [0, 0, 1, -E[2]],
  99. [0, 0, 0, 1]]
  100. ## end old
  101. return np.dot(Mr, Mt)
  102. def persp_transformation(zfront, zback):
  103. a = (zfront+zback)/(zfront-zback)
  104. b = -2*(zfront*zback)/(zfront-zback)
  105. return np.array([[1, 0, 0, 0],
  106. [0, 1, 0, 0],
  107. [0, 0, a, b],
  108. [0, 0, -1, 0]])
  109. def ortho_transformation(zfront, zback):
  110. # note: w component in the resulting vector will be (zback-zfront), not 1
  111. a = -(zfront + zback)
  112. b = -(zfront - zback)
  113. return np.array([[2, 0, 0, 0],
  114. [0, 2, 0, 0],
  115. [0, 0, -2, 0],
  116. [0, 0, a, b]])
  117. def _proj_transform_vec(vec, M):
  118. vecw = np.dot(M, vec)
  119. w = vecw[3]
  120. # clip here..
  121. txs, tys, tzs = vecw[0]/w, vecw[1]/w, vecw[2]/w
  122. return txs, tys, tzs
  123. @cbook.deprecated("3.1")
  124. def proj_transform_vec(vec, M):
  125. return _proj_transform_vec(vec, M)
  126. def _proj_transform_vec_clip(vec, M):
  127. vecw = np.dot(M, vec)
  128. w = vecw[3]
  129. # clip here.
  130. txs, tys, tzs = vecw[0] / w, vecw[1] / w, vecw[2] / w
  131. tis = (0 <= vecw[0]) & (vecw[0] <= 1) & (0 <= vecw[1]) & (vecw[1] <= 1)
  132. if np.any(tis):
  133. tis = vecw[1] < 1
  134. return txs, tys, tzs, tis
  135. @cbook.deprecated("3.1")
  136. def proj_transform_vec_clip(vec, M):
  137. return _proj_transform_vec_clip(vec, M)
  138. def inv_transform(xs, ys, zs, M):
  139. iM = linalg.inv(M)
  140. vec = _vec_pad_ones(xs, ys, zs)
  141. vecr = np.dot(iM, vec)
  142. try:
  143. vecr = vecr / vecr[3]
  144. except OverflowError:
  145. pass
  146. return vecr[0], vecr[1], vecr[2]
  147. def _vec_pad_ones(xs, ys, zs):
  148. return np.array([xs, ys, zs, np.ones_like(xs)])
  149. @cbook.deprecated("3.1")
  150. def vec_pad_ones(xs, ys, zs):
  151. return _vec_pad_ones(xs, ys, zs)
  152. def proj_transform(xs, ys, zs, M):
  153. """
  154. Transform the points by the projection matrix
  155. """
  156. vec = _vec_pad_ones(xs, ys, zs)
  157. return _proj_transform_vec(vec, M)
  158. transform = proj_transform
  159. def proj_transform_clip(xs, ys, zs, M):
  160. """
  161. Transform the points by the projection matrix
  162. and return the clipping result
  163. returns txs, tys, tzs, tis
  164. """
  165. vec = _vec_pad_ones(xs, ys, zs)
  166. return _proj_transform_vec_clip(vec, M)
  167. def proj_points(points, M):
  168. return np.column_stack(proj_trans_points(points, M))
  169. def proj_trans_points(points, M):
  170. xs, ys, zs = zip(*points)
  171. return proj_transform(xs, ys, zs, M)
  172. @cbook.deprecated("3.1")
  173. def proj_trans_clip_points(points, M):
  174. xs, ys, zs = zip(*points)
  175. return proj_transform_clip(xs, ys, zs, M)
  176. def rot_x(V, alpha):
  177. cosa, sina = np.cos(alpha), np.sin(alpha)
  178. M1 = np.array([[1, 0, 0, 0],
  179. [0, cosa, -sina, 0],
  180. [0, sina, cosa, 0],
  181. [0, 0, 0, 1]])
  182. return np.dot(M1, V)