Coverage for /usr/lib/python3/dist-packages/mpl_toolkits/mplot3d/proj3d.py: 65%

1""" 

2Various transforms used for by the 3D code 

3""" 

4 

5import numpy as np 

6import numpy.linalg as linalg 

7 

8 

9def _line2d_seg_dist(p1, p2, p0): 

10 """ 

11 Return the distance(s) from line defined by p1 - p2 to point(s) p0. 

12 

13 p0[0] = x(s) 

14 p0[1] = y(s) 

15 

16 intersection point p = p1 + u*(p2-p1) 

17 and intersection point lies within segment if u is between 0 and 1. 

18 

19 If p1 and p2 are identical, the distance between them and p0 is returned. 

20 """ 

21 

22 x01 = np.asarray(p0[0]) - p1[0] 

23 y01 = np.asarray(p0[1]) - p1[1] 

24 if np.all(p1[0:2] == p2[0:2]): 

25 return np.hypot(x01, y01) 

26 

27 x21 = p2[0] - p1[0] 

28 y21 = p2[1] - p1[1] 

29 u = (x01*x21 + y01*y21) / (x21**2 + y21**2) 

30 u = np.clip(u, 0, 1) 

31 d = np.hypot(x01 - u*x21, y01 - u*y21) 

32 

33 return d 

34 

35 

36def world_transformation(xmin, xmax, 

37 ymin, ymax, 

38 zmin, zmax, pb_aspect=None): 

39 """ 

40 Produce a matrix that scales homogeneous coords in the specified ranges 

41 to [0, 1], or [0, pb_aspect[i]] if the plotbox aspect ratio is specified. 

42 """ 

43 dx = xmax - xmin 

44 dy = ymax - ymin 

45 dz = zmax - zmin 

46 if pb_aspect is not None: 

47 ax, ay, az = pb_aspect 

48 dx /= ax 

49 dy /= ay 

50 dz /= az 

51 

52 return np.array([[1/dx, 0, 0, -xmin/dx], 

53 [0, 1/dy, 0, -ymin/dy], 

54 [0, 0, 1/dz, -zmin/dz], 

55 [0, 0, 0, 1]]) 

56 

57 

58def rotation_about_vector(v, angle): 

59 """ 

60 Produce a rotation matrix for an angle in radians about a vector. 

61 """ 

62 vx, vy, vz = v / np.linalg.norm(v) 

63 s = np.sin(angle) 

64 c = np.cos(angle) 

65 t = 2*np.sin(angle/2)**2 # more numerically stable than t = 1-c 

66 

67 R = np.array([ 

68 [t*vx*vx + c, t*vx*vy - vz*s, t*vx*vz + vy*s], 

69 [t*vy*vx + vz*s, t*vy*vy + c, t*vy*vz - vx*s], 

70 [t*vz*vx - vy*s, t*vz*vy + vx*s, t*vz*vz + c]]) 

71 

72 return R 

73 

74 

75def view_transformation(E, R, V, roll): 

76 n = (E - R) 

77 n = n/np.linalg.norm(n) 

78 u = np.cross(V, n) 

79 u = u/np.linalg.norm(u) 

80 v = np.cross(n, u) # Will be a unit vector 

81 

82 # Save some computation for the default roll=0 

83 if roll != 0: 

84 # A positive rotation of the camera is a negative rotation of the world 

85 Rroll = rotation_about_vector(n, -roll) 

86 u = np.dot(Rroll, u) 

87 v = np.dot(Rroll, v) 

88 

89 Mr = np.eye(4) 

90 Mt = np.eye(4) 

91 Mr[:3, :3] = [u, v, n] 

92 Mt[:3, -1] = -E 

93 

94 return np.dot(Mr, Mt) 

95 

96 

97def persp_transformation(zfront, zback, focal_length): 

98 e = focal_length 

99 a = 1 # aspect ratio 

100 b = (zfront+zback)/(zfront-zback) 

101 c = -2*(zfront*zback)/(zfront-zback) 

102 proj_matrix = np.array([[e, 0, 0, 0], 

103 [0, e/a, 0, 0], 

104 [0, 0, b, c], 

105 [0, 0, -1, 0]]) 

106 return proj_matrix 

107 

108 

109def 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 proj_matrix = np.array([[2, 0, 0, 0], 

114 [0, 2, 0, 0], 

115 [0, 0, -2, 0], 

116 [0, 0, a, b]]) 

117 return proj_matrix 

118 

119 

120def _proj_transform_vec(vec, M): 

121 vecw = np.dot(M, vec) 

122 w = vecw[3] 

123 # clip here.. 

124 txs, tys, tzs = vecw[0]/w, vecw[1]/w, vecw[2]/w 

125 return txs, tys, tzs 

126 

127 

128def _proj_transform_vec_clip(vec, M): 

129 vecw = np.dot(M, vec) 

130 w = vecw[3] 

131 # clip here. 

132 txs, tys, tzs = vecw[0] / w, vecw[1] / w, vecw[2] / w 

133 tis = (0 <= vecw[0]) & (vecw[0] <= 1) & (0 <= vecw[1]) & (vecw[1] <= 1) 

134 if np.any(tis): 

135 tis = vecw[1] < 1 

136 return txs, tys, tzs, tis 

137 

138 

139def inv_transform(xs, ys, zs, M): 

140 iM = linalg.inv(M) 

141 vec = _vec_pad_ones(xs, ys, zs) 

142 vecr = np.dot(iM, vec) 

143 try: 

144 vecr = vecr / vecr[3] 

145 except OverflowError: 

146 pass 

147 return vecr[0], vecr[1], vecr[2] 

148 

149 

150def _vec_pad_ones(xs, ys, zs): 

151 return np.array([xs, ys, zs, np.ones_like(xs)]) 

152 

153 

154def proj_transform(xs, ys, zs, M): 

155 """ 

156 Transform the points by the projection matrix 

157 """ 

158 vec = _vec_pad_ones(xs, ys, zs) 

159 return _proj_transform_vec(vec, M) 

160 

161 

162transform = proj_transform 

163 

164 

165def proj_transform_clip(xs, ys, zs, M): 

166 """ 

167 Transform the points by the projection matrix 

168 and return the clipping result 

169 returns txs, tys, tzs, tis 

170 """ 

171 vec = _vec_pad_ones(xs, ys, zs) 

172 return _proj_transform_vec_clip(vec, M) 

173 

174 

175def proj_points(points, M): 

176 return np.column_stack(proj_trans_points(points, M)) 

177 

178 

179def proj_trans_points(points, M): 

180 xs, ys, zs = zip(*points) 

181 return proj_transform(xs, ys, zs, M) 

182 

183 

184def rot_x(V, alpha): 

185 cosa, sina = np.cos(alpha), np.sin(alpha) 

186 M1 = np.array([[1, 0, 0, 0], 

187 [0, cosa, -sina, 0], 

188 [0, sina, cosa, 0], 

189 [0, 0, 0, 1]]) 

190 return np.dot(M1, V)