-
Notifications
You must be signed in to change notification settings - Fork 0
/
curvature.py
92 lines (75 loc) · 2.98 KB
/
curvature.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
import torch
from .diff import grad, div
# Compute the mean curvature of f at v (v: points on the surface f=0)
def meanCurvature(f, v):
if not torch.is_tensor(v):
# array.copy() is to prevent the case
# where the numpy array has a negative stride
vt = torch.tensor(v.copy())
vt.requires_grad = True
# eval the implicit at the points
ft = f(vt)
# derivatives
grad_ft = grad(ft, vt)
norm_grad_ft = torch.linalg.norm(grad_ft, 2, dim=1)
normalized_grad_ft = torch.zeros_like(grad_ft)
normalized_grad_ft[:, 0] = grad_ft[:, 0] / norm_grad_ft
normalized_grad_ft[:, 1] = grad_ft[:, 1] / norm_grad_ft
normalized_grad_ft[:, 2] = grad_ft[:, 2] / norm_grad_ft
mean_curvature = -0.5 * div(normalized_grad_ft, vt)
mc = mean_curvature[:, 0].detach().cpu().numpy()
return mc
def GaussianCurvature(f, v):
if not torch.is_tensor(v):
# array.copy() is to prevent the case
# where the numpy array has a negative stride
vt = torch.tensor(v.copy())
vt.requires_grad = True
# eval the implicit at the points
ft = f(vt)
# derivatives
grad_ft = grad(ft, vt)
norm_grad_ft = torch.linalg.norm(grad_ft, 2, dim=1)
norm_grad_ft4 = norm_grad_ft**4
grad_ft = grad_ft.reshape((grad_ft.shape[0], grad_ft.shape[1], 1))
# Hessian
H_ft = torch.zeros((v.shape[0], 3, 3))
tmp0 = grad(grad_ft[:, 0], vt)
H_ft[:, 0, 0] = tmp0[:, 0]
H_ft[:, 0, 1] = tmp0[:, 1]
H_ft[:, 0, 2] = tmp0[:, 2]
tmp1 = grad(grad_ft[:, 1], vt)
H_ft[:, 1, 0] = tmp1[:, 0]
H_ft[:, 1, 1] = tmp1[:, 1]
H_ft[:, 1, 2] = tmp1[:, 2]
tmp2 = grad(grad_ft[:, 2], vt)
H_ft[:, 2, 0] = tmp2[:, 0]
H_ft[:, 2, 1] = tmp2[:, 1]
H_ft[:, 2, 2] = tmp2[:, 2]
# adjoint of Hessian
Hstar_f = torch.zeros_like(H_ft)
Hstar_f[:, 0,
0] = H_ft[:, 1, 1] * H_ft[:, 2, 2] - H_ft[:, 1, 2] * H_ft[:, 2, 1]
Hstar_f[:, 1,
1] = H_ft[:, 0, 0] * H_ft[:, 2, 2] - H_ft[:, 0, 2] * H_ft[:, 2, 0]
Hstar_f[:, 2,
2] = H_ft[:, 0, 0] * H_ft[:, 1, 1] - H_ft[:, 0, 1] * H_ft[:, 1, 0]
Hstar_f[:, 1,
0] = H_ft[:, 0, 2] * H_ft[:, 2, 1] - H_ft[:, 0, 1] * H_ft[:, 2, 2]
Hstar_f[:, 2,
0] = H_ft[:, 0, 1] * H_ft[:, 1, 2] - H_ft[:, 0, 2] * H_ft[:, 1, 1]
Hstar_f[:, 0,
1] = H_ft[:, 1, 2] * H_ft[:, 2, 0] - H_ft[:, 1, 0] * H_ft[:, 2, 2]
Hstar_f[:, 2,
1] = H_ft[:, 1, 0] * H_ft[:, 0, 2] - H_ft[:, 0, 0] * H_ft[:, 1, 2]
Hstar_f[:, 0,
2] = H_ft[:, 1, 0] * H_ft[:, 2, 1] - H_ft[:, 1, 1] * H_ft[:, 2, 0]
Hstar_f[:, 1,
2] = H_ft[:, 0, 1] * H_ft[:, 2, 0] - H_ft[:, 0, 0] * H_ft[:, 2, 1]
tmp1 = torch.matmul(Hstar_f, grad_ft)
nv = grad_ft.shape[0]
tmp2 = torch.bmm(grad_ft.view(nv, 1, 3), tmp1.view(nv, 3, 1))
tmp3 = tmp2.reshape((nv))
tmp4 = tmp3 / norm_grad_ft4
Kg = tmp4.detach().cpu().numpy()
return Kg