proj3d.py 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175
  1. """
  2. Various transforms used for by the 3D code
  3. """
  4. import numpy as np
  5. import numpy.linalg as linalg
  6. def _line2d_seg_dist(p1, p2, p0):
  7. """
  8. Return the distance(s) from line defined by p1 - p2 to point(s) p0.
  9. p0[0] = x(s)
  10. p0[1] = y(s)
  11. intersection point p = p1 + u*(p2-p1)
  12. and intersection point lies within segment if u is between 0 and 1
  13. """
  14. x21 = p2[0] - p1[0]
  15. y21 = p2[1] - p1[1]
  16. x01 = np.asarray(p0[0]) - p1[0]
  17. y01 = np.asarray(p0[1]) - p1[1]
  18. u = (x01*x21 + y01*y21) / (x21**2 + y21**2)
  19. u = np.clip(u, 0, 1)
  20. d = np.hypot(x01 - u*x21, y01 - u*y21)
  21. return d
  22. def world_transformation(xmin, xmax,
  23. ymin, ymax,
  24. zmin, zmax, pb_aspect=None):
  25. """
  26. Produce a matrix that scales homogeneous coords in the specified ranges
  27. to [0, 1], or [0, pb_aspect[i]] if the plotbox aspect ratio is specified.
  28. """
  29. dx = xmax - xmin
  30. dy = ymax - ymin
  31. dz = zmax - zmin
  32. if pb_aspect is not None:
  33. ax, ay, az = pb_aspect
  34. dx /= ax
  35. dy /= ay
  36. dz /= az
  37. return np.array([[1/dx, 0, 0, -xmin/dx],
  38. [0, 1/dy, 0, -ymin/dy],
  39. [0, 0, 1/dz, -zmin/dz],
  40. [0, 0, 0, 1]])
  41. def view_transformation(E, R, V):
  42. n = (E - R)
  43. ## new
  44. # n /= np.linalg.norm(n)
  45. # u = np.cross(V, n)
  46. # u /= np.linalg.norm(u)
  47. # v = np.cross(n, u)
  48. # Mr = np.diag([1.] * 4)
  49. # Mt = np.diag([1.] * 4)
  50. # Mr[:3,:3] = u, v, n
  51. # Mt[:3,-1] = -E
  52. ## end new
  53. ## old
  54. n = n / np.linalg.norm(n)
  55. u = np.cross(V, n)
  56. u = u / np.linalg.norm(u)
  57. v = np.cross(n, u)
  58. Mr = [[u[0], u[1], u[2], 0],
  59. [v[0], v[1], v[2], 0],
  60. [n[0], n[1], n[2], 0],
  61. [0, 0, 0, 1]]
  62. #
  63. Mt = [[1, 0, 0, -E[0]],
  64. [0, 1, 0, -E[1]],
  65. [0, 0, 1, -E[2]],
  66. [0, 0, 0, 1]]
  67. ## end old
  68. return np.dot(Mr, Mt)
  69. def persp_transformation(zfront, zback):
  70. a = (zfront+zback)/(zfront-zback)
  71. b = -2*(zfront*zback)/(zfront-zback)
  72. return np.array([[1, 0, 0, 0],
  73. [0, 1, 0, 0],
  74. [0, 0, a, b],
  75. [0, 0, -1, 0]])
  76. def ortho_transformation(zfront, zback):
  77. # note: w component in the resulting vector will be (zback-zfront), not 1
  78. a = -(zfront + zback)
  79. b = -(zfront - zback)
  80. return np.array([[2, 0, 0, 0],
  81. [0, 2, 0, 0],
  82. [0, 0, -2, 0],
  83. [0, 0, a, b]])
  84. def _proj_transform_vec(vec, M):
  85. vecw = np.dot(M, vec)
  86. w = vecw[3]
  87. # clip here..
  88. txs, tys, tzs = vecw[0]/w, vecw[1]/w, vecw[2]/w
  89. return txs, tys, tzs
  90. def _proj_transform_vec_clip(vec, M):
  91. vecw = np.dot(M, vec)
  92. w = vecw[3]
  93. # clip here.
  94. txs, tys, tzs = vecw[0] / w, vecw[1] / w, vecw[2] / w
  95. tis = (0 <= vecw[0]) & (vecw[0] <= 1) & (0 <= vecw[1]) & (vecw[1] <= 1)
  96. if np.any(tis):
  97. tis = vecw[1] < 1
  98. return txs, tys, tzs, tis
  99. def inv_transform(xs, ys, zs, M):
  100. iM = linalg.inv(M)
  101. vec = _vec_pad_ones(xs, ys, zs)
  102. vecr = np.dot(iM, vec)
  103. try:
  104. vecr = vecr / vecr[3]
  105. except OverflowError:
  106. pass
  107. return vecr[0], vecr[1], vecr[2]
  108. def _vec_pad_ones(xs, ys, zs):
  109. return np.array([xs, ys, zs, np.ones_like(xs)])
  110. def proj_transform(xs, ys, zs, M):
  111. """
  112. Transform the points by the projection matrix
  113. """
  114. vec = _vec_pad_ones(xs, ys, zs)
  115. return _proj_transform_vec(vec, M)
  116. transform = proj_transform
  117. def proj_transform_clip(xs, ys, zs, M):
  118. """
  119. Transform the points by the projection matrix
  120. and return the clipping result
  121. returns txs, tys, tzs, tis
  122. """
  123. vec = _vec_pad_ones(xs, ys, zs)
  124. return _proj_transform_vec_clip(vec, M)
  125. def proj_points(points, M):
  126. return np.column_stack(proj_trans_points(points, M))
  127. def proj_trans_points(points, M):
  128. xs, ys, zs = zip(*points)
  129. return proj_transform(xs, ys, zs, M)
  130. def rot_x(V, alpha):
  131. cosa, sina = np.cos(alpha), np.sin(alpha)
  132. M1 = np.array([[1, 0, 0, 0],
  133. [0, cosa, -sina, 0],
  134. [0, sina, cosa, 0],
  135. [0, 0, 0, 1]])
  136. return np.dot(M1, V)