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
    # ------------------
    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