YoonaAI commited on
Commit
f7266a6
·
1 Parent(s): a2a65b6

Upload geometry.py

Browse files
Files changed (1) hide show
  1. lib/pymaf/utils/geometry.py +452 -0
lib/pymaf/utils/geometry.py ADDED
@@ -0,0 +1,452 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ import torch
2
+ import numpy as np
3
+ from torch.nn import functional as F
4
+ """
5
+ Useful geometric operations, e.g. Perspective projection and a differentiable Rodrigues formula
6
+ Parts of the code are taken from https://github.com/MandyMo/pytorch_HMR
7
+ """
8
+
9
+
10
+ def batch_rodrigues(theta):
11
+ """Convert axis-angle representation to rotation matrix.
12
+ Args:
13
+ theta: size = [B, 3]
14
+ Returns:
15
+ Rotation matrix corresponding to the quaternion -- size = [B, 3, 3]
16
+ """
17
+ l1norm = torch.norm(theta + 1e-8, p=2, dim=1)
18
+ angle = torch.unsqueeze(l1norm, -1)
19
+ normalized = torch.div(theta, angle)
20
+ angle = angle * 0.5
21
+ v_cos = torch.cos(angle)
22
+ v_sin = torch.sin(angle)
23
+ quat = torch.cat([v_cos, v_sin * normalized], dim=1)
24
+ return quat_to_rotmat(quat)
25
+
26
+
27
+ def quat_to_rotmat(quat):
28
+ """Convert quaternion coefficients to rotation matrix.
29
+ Args:
30
+ quat: size = [B, 4] 4 <===>(w, x, y, z)
31
+ Returns:
32
+ Rotation matrix corresponding to the quaternion -- size = [B, 3, 3]
33
+ """
34
+ norm_quat = quat
35
+ norm_quat = norm_quat / norm_quat.norm(p=2, dim=1, keepdim=True)
36
+ w, x, y, z = norm_quat[:, 0], norm_quat[:, 1], norm_quat[:,
37
+ 2], norm_quat[:,
38
+ 3]
39
+
40
+ B = quat.size(0)
41
+
42
+ w2, x2, y2, z2 = w.pow(2), x.pow(2), y.pow(2), z.pow(2)
43
+ wx, wy, wz = w * x, w * y, w * z
44
+ xy, xz, yz = x * y, x * z, y * z
45
+
46
+ rotMat = torch.stack([
47
+ w2 + x2 - y2 - z2, 2 * xy - 2 * wz, 2 * wy + 2 * xz, 2 * wz + 2 * xy,
48
+ w2 - x2 + y2 - z2, 2 * yz - 2 * wx, 2 * xz - 2 * wy, 2 * wx + 2 * yz,
49
+ w2 - x2 - y2 + z2
50
+ ],
51
+ dim=1).view(B, 3, 3)
52
+ return rotMat
53
+
54
+
55
+ def rotation_matrix_to_angle_axis(rotation_matrix):
56
+ """
57
+ This function is borrowed from https://github.com/kornia/kornia
58
+
59
+ Convert 3x4 rotation matrix to Rodrigues vector
60
+
61
+ Args:
62
+ rotation_matrix (Tensor): rotation matrix.
63
+
64
+ Returns:
65
+ Tensor: Rodrigues vector transformation.
66
+
67
+ Shape:
68
+ - Input: :math:`(N, 3, 4)`
69
+ - Output: :math:`(N, 3)`
70
+
71
+ Example:
72
+ >>> input = torch.rand(2, 3, 4) # Nx4x4
73
+ >>> output = tgm.rotation_matrix_to_angle_axis(input) # Nx3
74
+ """
75
+ if rotation_matrix.shape[1:] == (3, 3):
76
+ rot_mat = rotation_matrix.reshape(-1, 3, 3)
77
+ hom = torch.tensor([0, 0, 1],
78
+ dtype=torch.float32,
79
+ device=rotation_matrix.device).reshape(
80
+ 1, 3, 1).expand(rot_mat.shape[0], -1, -1)
81
+ rotation_matrix = torch.cat([rot_mat, hom], dim=-1)
82
+
83
+ quaternion = rotation_matrix_to_quaternion(rotation_matrix)
84
+ aa = quaternion_to_angle_axis(quaternion)
85
+ aa[torch.isnan(aa)] = 0.0
86
+ return aa
87
+
88
+
89
+ def quaternion_to_angle_axis(quaternion: torch.Tensor) -> torch.Tensor:
90
+ """
91
+ This function is borrowed from https://github.com/kornia/kornia
92
+
93
+ Convert quaternion vector to angle axis of rotation.
94
+
95
+ Adapted from ceres C++ library: ceres-solver/include/ceres/rotation.h
96
+
97
+ Args:
98
+ quaternion (torch.Tensor): tensor with quaternions.
99
+
100
+ Return:
101
+ torch.Tensor: tensor with angle axis of rotation.
102
+
103
+ Shape:
104
+ - Input: :math:`(*, 4)` where `*` means, any number of dimensions
105
+ - Output: :math:`(*, 3)`
106
+
107
+ Example:
108
+ >>> quaternion = torch.rand(2, 4) # Nx4
109
+ >>> angle_axis = tgm.quaternion_to_angle_axis(quaternion) # Nx3
110
+ """
111
+ if not torch.is_tensor(quaternion):
112
+ raise TypeError("Input type is not a torch.Tensor. Got {}".format(
113
+ type(quaternion)))
114
+
115
+ if not quaternion.shape[-1] == 4:
116
+ raise ValueError(
117
+ "Input must be a tensor of shape Nx4 or 4. Got {}".format(
118
+ quaternion.shape))
119
+ # unpack input and compute conversion
120
+ q1: torch.Tensor = quaternion[..., 1]
121
+ q2: torch.Tensor = quaternion[..., 2]
122
+ q3: torch.Tensor = quaternion[..., 3]
123
+ sin_squared_theta: torch.Tensor = q1 * q1 + q2 * q2 + q3 * q3
124
+
125
+ sin_theta: torch.Tensor = torch.sqrt(sin_squared_theta)
126
+ cos_theta: torch.Tensor = quaternion[..., 0]
127
+ two_theta: torch.Tensor = 2.0 * torch.where(
128
+ cos_theta < 0.0, torch.atan2(-sin_theta, -cos_theta),
129
+ torch.atan2(sin_theta, cos_theta))
130
+
131
+ k_pos: torch.Tensor = two_theta / sin_theta
132
+ k_neg: torch.Tensor = 2.0 * torch.ones_like(sin_theta)
133
+ k: torch.Tensor = torch.where(sin_squared_theta > 0.0, k_pos, k_neg)
134
+
135
+ angle_axis: torch.Tensor = torch.zeros_like(quaternion)[..., :3]
136
+ angle_axis[..., 0] += q1 * k
137
+ angle_axis[..., 1] += q2 * k
138
+ angle_axis[..., 2] += q3 * k
139
+ return angle_axis
140
+
141
+
142
+ def rotation_matrix_to_quaternion(rotation_matrix, eps=1e-6):
143
+ """
144
+ This function is borrowed from https://github.com/kornia/kornia
145
+
146
+ Convert 3x4 rotation matrix to 4d quaternion vector
147
+
148
+ This algorithm is based on algorithm described in
149
+ https://github.com/KieranWynn/pyquaternion/blob/master/pyquaternion/quaternion.py#L201
150
+
151
+ Args:
152
+ rotation_matrix (Tensor): the rotation matrix to convert.
153
+
154
+ Return:
155
+ Tensor: the rotation in quaternion
156
+
157
+ Shape:
158
+ - Input: :math:`(N, 3, 4)`
159
+ - Output: :math:`(N, 4)`
160
+
161
+ Example:
162
+ >>> input = torch.rand(4, 3, 4) # Nx3x4
163
+ >>> output = tgm.rotation_matrix_to_quaternion(input) # Nx4
164
+ """
165
+ if not torch.is_tensor(rotation_matrix):
166
+ raise TypeError("Input type is not a torch.Tensor. Got {}".format(
167
+ type(rotation_matrix)))
168
+
169
+ if len(rotation_matrix.shape) > 3:
170
+ raise ValueError(
171
+ "Input size must be a three dimensional tensor. Got {}".format(
172
+ rotation_matrix.shape))
173
+ if not rotation_matrix.shape[-2:] == (3, 4):
174
+ raise ValueError(
175
+ "Input size must be a N x 3 x 4 tensor. Got {}".format(
176
+ rotation_matrix.shape))
177
+
178
+ rmat_t = torch.transpose(rotation_matrix, 1, 2)
179
+
180
+ mask_d2 = rmat_t[:, 2, 2] < eps
181
+
182
+ mask_d0_d1 = rmat_t[:, 0, 0] > rmat_t[:, 1, 1]
183
+ mask_d0_nd1 = rmat_t[:, 0, 0] < -rmat_t[:, 1, 1]
184
+
185
+ t0 = 1 + rmat_t[:, 0, 0] - rmat_t[:, 1, 1] - rmat_t[:, 2, 2]
186
+ q0 = torch.stack([
187
+ rmat_t[:, 1, 2] - rmat_t[:, 2, 1], t0,
188
+ rmat_t[:, 0, 1] + rmat_t[:, 1, 0], rmat_t[:, 2, 0] + rmat_t[:, 0, 2]
189
+ ], -1)
190
+ t0_rep = t0.repeat(4, 1).t()
191
+
192
+ t1 = 1 - rmat_t[:, 0, 0] + rmat_t[:, 1, 1] - rmat_t[:, 2, 2]
193
+ q1 = torch.stack([
194
+ rmat_t[:, 2, 0] - rmat_t[:, 0, 2], rmat_t[:, 0, 1] + rmat_t[:, 1, 0],
195
+ t1, rmat_t[:, 1, 2] + rmat_t[:, 2, 1]
196
+ ], -1)
197
+ t1_rep = t1.repeat(4, 1).t()
198
+
199
+ t2 = 1 - rmat_t[:, 0, 0] - rmat_t[:, 1, 1] + rmat_t[:, 2, 2]
200
+ q2 = torch.stack([
201
+ rmat_t[:, 0, 1] - rmat_t[:, 1, 0], rmat_t[:, 2, 0] + rmat_t[:, 0, 2],
202
+ rmat_t[:, 1, 2] + rmat_t[:, 2, 1], t2
203
+ ], -1)
204
+ t2_rep = t2.repeat(4, 1).t()
205
+
206
+ t3 = 1 + rmat_t[:, 0, 0] + rmat_t[:, 1, 1] + rmat_t[:, 2, 2]
207
+ q3 = torch.stack([
208
+ t3, rmat_t[:, 1, 2] - rmat_t[:, 2, 1],
209
+ rmat_t[:, 2, 0] - rmat_t[:, 0, 2], rmat_t[:, 0, 1] - rmat_t[:, 1, 0]
210
+ ], -1)
211
+ t3_rep = t3.repeat(4, 1).t()
212
+
213
+ mask_c0 = mask_d2 * mask_d0_d1
214
+ mask_c1 = mask_d2 * ~mask_d0_d1
215
+ mask_c2 = ~mask_d2 * mask_d0_nd1
216
+ mask_c3 = ~mask_d2 * ~mask_d0_nd1
217
+ mask_c0 = mask_c0.view(-1, 1).type_as(q0)
218
+ mask_c1 = mask_c1.view(-1, 1).type_as(q1)
219
+ mask_c2 = mask_c2.view(-1, 1).type_as(q2)
220
+ mask_c3 = mask_c3.view(-1, 1).type_as(q3)
221
+
222
+ q = q0 * mask_c0 + q1 * mask_c1 + q2 * mask_c2 + q3 * mask_c3
223
+ q /= torch.sqrt(t0_rep * mask_c0 + t1_rep * mask_c1 + # noqa
224
+ t2_rep * mask_c2 + t3_rep * mask_c3) # noqa
225
+ q *= 0.5
226
+ return q
227
+
228
+
229
+ def rot6d_to_rotmat(x):
230
+ """Convert 6D rotation representation to 3x3 rotation matrix.
231
+ Based on Zhou et al., "On the Continuity of Rotation Representations in Neural Networks", CVPR 2019
232
+ Input:
233
+ (B,6) Batch of 6-D rotation representations
234
+ Output:
235
+ (B,3,3) Batch of corresponding rotation matrices
236
+ """
237
+ x = x.view(-1, 3, 2)
238
+ a1 = x[:, :, 0]
239
+ a2 = x[:, :, 1]
240
+ b1 = F.normalize(a1)
241
+ b2 = F.normalize(a2 - torch.einsum('bi,bi->b', b1, a2).unsqueeze(-1) * b1)
242
+ b3 = torch.cross(b1, b2)
243
+ return torch.stack((b1, b2, b3), dim=-1)
244
+
245
+
246
+ def projection(pred_joints, pred_camera, retain_z=False):
247
+ pred_cam_t = torch.stack([
248
+ pred_camera[:, 1], pred_camera[:, 2], 2 * 5000. /
249
+ (224. * pred_camera[:, 0] + 1e-9)
250
+ ],
251
+ dim=-1)
252
+ batch_size = pred_joints.shape[0]
253
+ camera_center = torch.zeros(batch_size, 2)
254
+ pred_keypoints_2d = perspective_projection(
255
+ pred_joints,
256
+ rotation=torch.eye(3).unsqueeze(0).expand(batch_size, -1,
257
+ -1).to(pred_joints.device),
258
+ translation=pred_cam_t,
259
+ focal_length=5000.,
260
+ camera_center=camera_center,
261
+ retain_z=retain_z)
262
+ # Normalize keypoints to [-1,1]
263
+ pred_keypoints_2d = pred_keypoints_2d / (224. / 2.)
264
+ return pred_keypoints_2d
265
+
266
+
267
+ def perspective_projection(points,
268
+ rotation,
269
+ translation,
270
+ focal_length,
271
+ camera_center,
272
+ retain_z=False):
273
+ """
274
+ This function computes the perspective projection of a set of points.
275
+ Input:
276
+ points (bs, N, 3): 3D points
277
+ rotation (bs, 3, 3): Camera rotation
278
+ translation (bs, 3): Camera translation
279
+ focal_length (bs,) or scalar: Focal length
280
+ camera_center (bs, 2): Camera center
281
+ """
282
+ batch_size = points.shape[0]
283
+ K = torch.zeros([batch_size, 3, 3], device=points.device)
284
+ K[:, 0, 0] = focal_length
285
+ K[:, 1, 1] = focal_length
286
+ K[:, 2, 2] = 1.
287
+ K[:, :-1, -1] = camera_center
288
+
289
+ # Transform points
290
+ points = torch.einsum('bij,bkj->bki', rotation, points)
291
+ points = points + translation.unsqueeze(1)
292
+
293
+ # Apply perspective distortion
294
+ projected_points = points / points[:, :, -1].unsqueeze(-1)
295
+
296
+ # Apply camera intrinsics
297
+ projected_points = torch.einsum('bij,bkj->bki', K, projected_points)
298
+
299
+ if retain_z:
300
+ return projected_points
301
+ else:
302
+ return projected_points[:, :, :-1]
303
+
304
+
305
+ def estimate_translation_np(S,
306
+ joints_2d,
307
+ joints_conf,
308
+ focal_length=5000,
309
+ img_size=224):
310
+ """Find camera translation that brings 3D joints S closest to 2D the corresponding joints_2d.
311
+ Input:
312
+ S: (25, 3) 3D joint locations
313
+ joints: (25, 3) 2D joint locations and confidence
314
+ Returns:
315
+ (3,) camera translation vector
316
+ """
317
+
318
+ num_joints = S.shape[0]
319
+ # focal length
320
+ f = np.array([focal_length, focal_length])
321
+ # optical center
322
+ center = np.array([img_size / 2., img_size / 2.])
323
+
324
+ # transformations
325
+ Z = np.reshape(np.tile(S[:, 2], (2, 1)).T, -1)
326
+ XY = np.reshape(S[:, 0:2], -1)
327
+ O = np.tile(center, num_joints)
328
+ F = np.tile(f, num_joints)
329
+ weight2 = np.reshape(np.tile(np.sqrt(joints_conf), (2, 1)).T, -1)
330
+
331
+ # least squares
332
+ Q = np.array([
333
+ F * np.tile(np.array([1, 0]), num_joints),
334
+ F * np.tile(np.array([0, 1]), num_joints),
335
+ O - np.reshape(joints_2d, -1)
336
+ ]).T
337
+ c = (np.reshape(joints_2d, -1) - O) * Z - F * XY
338
+
339
+ # weighted least squares
340
+ W = np.diagflat(weight2)
341
+ Q = np.dot(W, Q)
342
+ c = np.dot(W, c)
343
+
344
+ # square matrix
345
+ A = np.dot(Q.T, Q)
346
+ b = np.dot(Q.T, c)
347
+
348
+ # solution
349
+ trans = np.linalg.solve(A, b)
350
+
351
+ return trans
352
+
353
+
354
+ def estimate_translation(S, joints_2d, focal_length=5000., img_size=224.):
355
+ """Find camera translation that brings 3D joints S closest to 2D the corresponding joints_2d.
356
+ Input:
357
+ S: (B, 49, 3) 3D joint locations
358
+ joints: (B, 49, 3) 2D joint locations and confidence
359
+ Returns:
360
+ (B, 3) camera translation vectors
361
+ """
362
+
363
+ device = S.device
364
+ # Use only joints 25:49 (GT joints)
365
+ S = S[:, 25:, :].cpu().numpy()
366
+ joints_2d = joints_2d[:, 25:, :].cpu().numpy()
367
+ joints_conf = joints_2d[:, :, -1]
368
+ joints_2d = joints_2d[:, :, :-1]
369
+ trans = np.zeros((S.shape[0], 3), dtype=np.float32)
370
+ # Find the translation for each example in the batch
371
+ for i in range(S.shape[0]):
372
+ S_i = S[i]
373
+ joints_i = joints_2d[i]
374
+ conf_i = joints_conf[i]
375
+ trans[i] = estimate_translation_np(S_i,
376
+ joints_i,
377
+ conf_i,
378
+ focal_length=focal_length,
379
+ img_size=img_size)
380
+ return torch.from_numpy(trans).to(device)
381
+
382
+
383
+ def Rot_y(angle, category='torch', prepend_dim=True, device=None):
384
+ '''Rotate around y-axis by angle
385
+ Args:
386
+ category: 'torch' or 'numpy'
387
+ prepend_dim: prepend an extra dimension
388
+ Return: Rotation matrix with shape [1, 3, 3] (prepend_dim=True)
389
+ '''
390
+ m = np.array([[np.cos(angle), 0., np.sin(angle)], [0., 1., 0.],
391
+ [-np.sin(angle), 0., np.cos(angle)]])
392
+ if category == 'torch':
393
+ if prepend_dim:
394
+ return torch.tensor(m, dtype=torch.float,
395
+ device=device).unsqueeze(0)
396
+ else:
397
+ return torch.tensor(m, dtype=torch.float, device=device)
398
+ elif category == 'numpy':
399
+ if prepend_dim:
400
+ return np.expand_dims(m, 0)
401
+ else:
402
+ return m
403
+ else:
404
+ raise ValueError("category must be 'torch' or 'numpy'")
405
+
406
+
407
+ def Rot_x(angle, category='torch', prepend_dim=True, device=None):
408
+ '''Rotate around x-axis by angle
409
+ Args:
410
+ category: 'torch' or 'numpy'
411
+ prepend_dim: prepend an extra dimension
412
+ Return: Rotation matrix with shape [1, 3, 3] (prepend_dim=True)
413
+ '''
414
+ m = np.array([[1., 0., 0.], [0., np.cos(angle), -np.sin(angle)],
415
+ [0., np.sin(angle), np.cos(angle)]])
416
+ if category == 'torch':
417
+ if prepend_dim:
418
+ return torch.tensor(m, dtype=torch.float,
419
+ device=device).unsqueeze(0)
420
+ else:
421
+ return torch.tensor(m, dtype=torch.float, device=device)
422
+ elif category == 'numpy':
423
+ if prepend_dim:
424
+ return np.expand_dims(m, 0)
425
+ else:
426
+ return m
427
+ else:
428
+ raise ValueError("category must be 'torch' or 'numpy'")
429
+
430
+
431
+ def Rot_z(angle, category='torch', prepend_dim=True, device=None):
432
+ '''Rotate around z-axis by angle
433
+ Args:
434
+ category: 'torch' or 'numpy'
435
+ prepend_dim: prepend an extra dimension
436
+ Return: Rotation matrix with shape [1, 3, 3] (prepend_dim=True)
437
+ '''
438
+ m = np.array([[np.cos(angle), -np.sin(angle), 0.],
439
+ [np.sin(angle), np.cos(angle), 0.], [0., 0., 1.]])
440
+ if category == 'torch':
441
+ if prepend_dim:
442
+ return torch.tensor(m, dtype=torch.float,
443
+ device=device).unsqueeze(0)
444
+ else:
445
+ return torch.tensor(m, dtype=torch.float, device=device)
446
+ elif category == 'numpy':
447
+ if prepend_dim:
448
+ return np.expand_dims(m, 0)
449
+ else:
450
+ return m
451
+ else:
452
+ raise ValueError("category must be 'torch' or 'numpy'")