# Untitled

unknown
plain_text
19 days ago
1.5 kB
3
Indexable
Never
```def shape_operator(m,v):
pos = m.positions()

# Insert code here which computes the shape operator matrix S that can be used to compute
# the principal curvatures kmin and kmax and the mean curvature H

# ------------------
for v in m.vertices():
normal = m.vertex_normal(v)
n1, n2, n3 = normal
if n1 == 0:
a_ = [1, n2, n3]
else:
a_ = [0, n2, n3]
# a_ = [-n1, n2, n3]
length= np.linalg.norm(a_-normal*(np.dot(a_,normal)))
if length == 0:
print("OH NO")
length = 1e-10
a = a_-normal*(np.dot(a_,normal))/length
b = np.cross(normal, a)
T = np.array((a,b)).T
A = np.zeros((len(m.circulate_vertex(v, 'v')), 3))
B = np.zeros(len(m.circulate_vertex(v, 'v')))
pi = pos[v]
for idx,pj_idx in enumerate(m.circulate_vertex(v, 'v')):
pj = pos[pj_idx]
uj, vj = T.T@(pj-pi)
h = normal@(pj-pi)
A[idx] = np.array(((uj**2)/2, uj*vj, (vj**2)/2))
B[idx] = h
# Solve the least squares problem
a, b, c = np.linalg.lstsq(A,B, rcond=None)[0]
# print(a,b,c)
S = -np.array([[a, b], [b, c]])
# print(S)
# print(np.linalg.eig(S))
dir1, dir2= np.linalg.eig(S)[0]
if dir1>dir2:
kmax = dir1
kmin = dir2
else:
kmax = dir2
kmin = dir1
H = (kmax+kmin)/2

return kmin, kmax, H```