property_topomesh_optimization.py 文件源码

def property_topomesh_cotangent_laplacian_smoothing_force(topomesh):

    compute_topomesh_property(topomesh,'vertices',degree=2)
    compute_topomesh_property(topomesh,'length',degree=1)

    triangle_vertices = topomesh.wisp_property('vertices',degree=2).values(list(topomesh.wisps(2)))
    rotated_triangle_vertices = np.transpose([triangle_vertices[:,2],triangle_vertices[:,0],triangle_vertices[:,1]])
    antirotated_triangle_vertices = np.transpose([triangle_vertices[:,1],triangle_vertices[:,2],triangle_vertices[:,0]])
    triangle_vertices = np.append(np.append(triangle_vertices,rotated_triangle_vertices,axis=0),antirotated_triangle_vertices,axis=0)

    edge_index_list = np.array([[1, 2],[0, 2],[0, 1]])
    triangle_edge_vertices = triangle_vertices[:,edge_index_list]

    triangle_edge_vectors = topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,1]) - topomesh.wisp_property('barycenter',degree=0).values(triangle_edge_vertices[...,0])
    #triangle_edge_lengths = np.power(np.sum(np.power(triangle_edge_vectors,2.0),axis=2),0.5)
    triangle_edge_lengths = np.linalg.norm(triangle_edge_vectors,axis=2)
    triangle_edge_directions = triangle_edge_vectors/triangle_edge_lengths[...,np.newaxis]

    triangle_perimeters = np.sum(triangle_edge_lengths,axis=1)
    triangle_areas = np.sqrt((triangle_perimeters/2.0)*(triangle_perimeters/2.0-triangle_edge_lengths[:,0])*(triangle_perimeters/2.0-triangle_edge_lengths[:,1])*(triangle_perimeters/2.0-triangle_edge_lengths[:,2]))

    triangle_cosines = np.zeros_like(triangle_edge_lengths,np.float32)
    triangle_cosines[:,0] = (triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2)/(2.0*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2])
    triangle_cosines[:,1] = (triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2)/(2.0*triangle_edge_lengths[:,2]*triangle_edge_lengths[:,0])
    triangle_cosines[:,2] = (triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2)/(2.0*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1])
    triangle_cosines[np.where(np.abs(triangle_cosines) < np.power(10.,-6))] = np.power(10.,-6)
    triangle_angles = np.arccos(triangle_cosines)

    triangle_sinuses = np.zeros_like(triangle_edge_lengths,np.float32)
    triangle_sinuses[:,0] = np.sqrt(np.array(1.0 - np.power(triangle_edge_lengths[:,1]**2+triangle_edge_lengths[:,2]**2-triangle_edge_lengths[:,0]**2,2.0)/np.power(2.0*triangle_edge_lengths[:,1]*triangle_edge_lengths[:,2],2.0),np.float16))
    triangle_sinuses[:,1] = np.sqrt(np.array(1.0 - np.power(triangle_edge_lengths[:,2]**2+triangle_edge_lengths[:,0]**2-triangle_edge_lengths[:,1]**2,2.0)/np.power(2.0*triangle_edge_lengths[:,2]*triangle_edge_lengths[:,0],2.0),np.float16))
    triangle_sinuses[:,2] = np.sqrt(np.array(1.0 - np.power(triangle_edge_lengths[:,0]**2+triangle_edge_lengths[:,1]**2-triangle_edge_lengths[:,2]**2,2.0)/np.power(2.0*triangle_edge_lengths[:,0]*triangle_edge_lengths[:,1],2.0),np.float16))
    triangle_sinuses[np.where(triangle_sinuses < np.power(10.,-6))] = np.power(10.,-6)

    triangle_cotangent_weights = (triangle_cosines/triangle_sinuses)
    #triangle_cotangent_weights = 1./triangle_edge_lengths
    #triangle_cotangent_weights = 0.5*(triangle_cosines/triangle_sinuses)/(triangle_edge_lengths + np.power(10.,-10))
    #triangle_cotangent_vectors = (triangle_cosines/triangle_sinuses)[...,np.newaxis] * triangle_edge_vectors/triangle_edge_lengths[:,np.newaxis]
    triangle_cotangent_vectors = triangle_cotangent_weights[...,np.newaxis] * triangle_edge_vectors
    #triangle_cotangent_vectors = triangle_cotangent_weights[...,np.newaxis] * triangle_edge_directions
    # triangle_cotangent_vectors = 1./np.tan(triangle_angles)[...,np.newaxis] * triangle_edge_vectors

    vertex_cotangent_force = np.transpose([nd.sum(triangle_cotangent_vectors[:,1,d]+triangle_cotangent_vectors[:,2,d],triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0)))) for d in xrange(3)])
    vertex_cotangent_sum = nd.sum(triangle_cotangent_weights[:,1] + triangle_cotangent_weights[:,2],triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0))))

    vertex_area = nd.sum(triangle_areas,triangle_vertices[:,0],index=np.array(list(topomesh.wisps(0))))

    #laplacian_force = vertex_cotangent_force
    laplacian_force = vertex_cotangent_force/vertex_cotangent_sum[...,np.newaxis]
    #laplacian_force = vertex_cotangent_force/(4.*vertex_area[...,np.newaxis])
    laplacian_force[np.where(np.isnan(laplacian_force))] = 0.0
    laplacian_force[np.where(np.isinf(laplacian_force))] = 0.0

    return laplacian_force