# Untitled

unknown
python
12 days ago
2.4 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

# ------------------
N = m.vertex_normal(v)
F = []
U = []
a_tilde = N
a_tilde[1] *= -1
projectionOnA = a_tilde - N*(np.dot(a_tilde,N))
a = (projectionOnA)/np.linalg.norm(projectionOnA)
b = np.cross(N,a)
pi = pos[v]

T = np.array([a,b])
neighbours = list(m.circulate_vertex(v,"v"))

U = np.zeros((len(neighbours), 3))
F = np.zeros((len(neighbours), 1))

for j, neighbor in enumerate(neighbours):
pj = pos[neighbor]
uv = T @ (pj - pi)
hj = np.dot(N, (pj - pi))
U[j, :] = [uv[0]**2, uv[0]*uv[1], uv[1]**2]
F[j] = hj

# Solve the least squares problem
coeffs, residuals, rank, s = np.linalg.lstsq(U, F, rcond=None)

#Get 2D Shape operator
S = -1*np.array([[coeffs[0][0], coeffs[1][0]], [coeffs[1][0], coeffs[2][0]]])

#Convert to 3D
S = T.T@S@T

eigVal,eigVec = np.linalg.eig(S)

#The three vectors
#1. Eigenvector closest to 0 -> normal
#2. Lowest -> min
#3. Highest -> max

#First we get the normal vector
normalIdx = np.argmin(np.abs(eigVal))
normalVector = eigVec[normalIdx]
normalValue = eigVal[normalIdx]

#Then we check if any of the other vectors are as close to 0
numericalError = max(1**(-1000),normalValue)
zeroValueIdx = []

for x in range(len(eigVal)):
absEigenValue = abs(eigVal[x])
if(absEigenValue <= numericalError):
zeroValueIdx += [x]

nrZeroValues = len(zeroValueIdx)
print(zeroValueIdx)

#Error
if(nrZeroValues == 0):
return 0/0
if(nrZeroValues == 1):
np.delete(eigVal,normalIdx)
np.delete(eigVec,normalIdx)
kmax = max(eigVal)
kmin = min(eigVal)
if(nrZeroValues == 2):
print("Two 0's - but we're just gonna ignore it instead of figuring out the principal curvature :D")
print(eigVal)
np.delete(eigVal,normalIdx)
np.delete(eigVec,normalIdx)
kmax = max(eigVal)
kmin = min(eigVal)

# Mean curvature
H = np.linalg.det(S)

return kmin, kmax, H```