Untitled

def normalize(v):
    return v / linalg.norm(v)


# The function below is a helper function. The idea is that you give it a mesh and
# three vertices from a triangle. i and j define the edge, and k is the vertex at which 
# you compute cotan of the angle. 
def cotan_angle(m,i,j,k):
    # Compute and return the cotan
    pos = m.positions()
    v1 = pos[i] - pos[k]
    v2 = pos[j] - pos[k]
    cot = np.dot(v1,v2)/np.linalg.norm(np.cross(v1,v2))
    return cot

# This function is passed a mesh and returns the cotan LBO 
def create_laplacian(m, symmetric=False):
    # symmetric is a flag that you can choose in order to make the matrix L symmetric or not
    no_verts = m.no_allocated_vertices()
    L = zeros((no_verts,no_verts)) # The cot Laplacian that we end up outputting 
    A = zeros((no_verts)) # Use this vector to store 1-ring areas

    for v in m.vertices():
        ring_points = list(m.circulate_vertex(v, 'v'))
        no_neighbours = (len(ring_points))
        ring_points.append(ring_points[0])
        ring_points.append(ring_points[1])
        sum = 0
        Ai = 0
        for f in m.circulate_vertex(v,"f"):
            Ai += m.area(f)
        A[v] = Ai
        Li = 0
        #First we get the cotan weights
        for n in range(1, no_neighbours+1):
            j = ring_points[n]

            cot_alpha = cotan_angle(m, v, ring_points[n],ring_points[n-1])
            cot_beta = cotan_angle(m, v, ring_points[n],ring_points[n+1])
            Lij = (cot_alpha+ cot_beta)/(2*Ai)
            L[v][j] = Lij
            Li += Lij

        L[v][v] = -Li


    
    
    return L

# Load the mesh and make safe copy of vertex positions
m = hmesh.load("bunnygtest.obj")
pos = m.positions()
pos_orig = array(pos)
N = m.no_allocated_vertices()

# Explicit smoothing.
#
# Explicit smoothing is simply a matrix multiplication. Try several step sizes and 
# note how the effect changes. Also, try to perform several iterations of smoothing.
#

# Your code here
# ---------------------------

alpha = 0.5
delta = 0.5*10**(-3)
t = 1
for x in range(t):
    L = create_laplacian(m,False)
    pos[:] = (pos + alpha*delta*L@pos)

jd.display(m,smooth=False)

# Resetting to noisy positions
pos[:] = pos_noisy