Untitled

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
    
    # Your code here
    # ------------------
    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