CG——Laplacian Surface Editing

参考文献：https://igl.ethz.ch/projects/Laplacian-mesh-processing/Laplacian-mesh-editing/laplacian-mesh-editing.pdf

https://people.eecs.berkeley.edu/~jrs/meshpapers/Sorkine.pdf

原理：这篇博客写的很好，通俗易懂：https://blog.csdn.net/qq_31804159/article/details/103409465

代码：

import numpy as np
import networkx as nx
import trimesh
import scipy.sparse

#function collection
def mesh_to_graph(mesh):
    edges = mesh.edges_unique
    verts = mesh.vertices
    g = nx.Graph()
    for i, pos in enumerate(verts):
        g.add_node(i, pos=pos)
    for edge in edges:
        g.add_edge(*edge)
    return g

def get_boundary(g, ctrl_points):
    boundary = []
    for i in range(len(ctrl_points)):
        p1 = ctrl_points[i]
        p2 = ctrl_points[(i+1)%len(ctrl_points)]
        path = nx.shortest_path(g, p1, p2)
        boundary += path[0:-1]
    return boundary

def rw_laplacian_matrix(g):
    lap = nx.laplacian_matrix(g)
    adj = nx.adjacency_matrix(g)
    inv_deg = scipy.sparse.diags(1/adj.dot(np.ones([adj.shape[0]])))
    return inv_deg.dot(lap)

def get_editable_vertices(g, boundary, manip_handle):
    g_part = g.copy()
    g_part.remove_nodes_from(boundary)
    return list(nx.node_connected_component(g_part, manip_handle))

def get_mesh_scene(mesh, boundary=None, editable_vertices=None):
    scene = [mesh]
#     if boundary:
#         scene.append(trimesh.load_path(mesh.vertices[boundary+[boundary[0]]]))
    if editable_vertices:
        scene.append(trimesh.points.PointCloud(mesh.vertices[editable_vertices + boundary]))
    return trimesh.Scene(scene)

#load mesh

mesh = trimesh.load('bunny.ply', process=False)
g = mesh_to_graph(mesh)
handles = {
    24092 : [-0.01134  ,  0.151374 , -0.0242688]
}
boundary_ctrl_points = [15617, 24120, 30216, 11236, 6973]
boundary = get_boundary(g, boundary_ctrl_points)
editable_verts = get_editable_vertices(g, boundary, list(handles.keys())[0])
trimesh.points.PointCloud(mesh.vertices[editable_verts + boundary])
get_mesh_scene(mesh, boundary, editable_verts).show(viewer='gl',smooth=False)

# laplacian
subgraph = g.subgraph(boundary + editable_verts)
old_id = list(subgraph)
grp_idx = {}
for node in subgraph.nodes:
    grp_idx[node] = len(grp_idx)
g2 = list(subgraph.nodes)
L = rw_laplacian_matrix(subgraph).todense()   #最初的L
V = np.matrix([subgraph.nodes[i]['pos'] for i in subgraph.nodes])
Delta = L.dot(V)
boundary_idx = [grp_idx[i] for i in boundary_ctrl_points]
for x in handles.keys():
    t = grp_idx[int(x)]
    boundary_idx.append(t)
N = L.shape[0]       
H = np.zeros((len(boundary_idx), N))
h = [] 
for i in range(len(boundary_idx)):
    idx = boundary_idx[i]
    H[i][idx] = 1
for x in boundary_ctrl_points:
    h.append( np.array(subgraph.nodes[x]['pos']) )
for x in handles.values():
    h.append( x )
h = np.matrix(h)
newL = np.vstack((L, H))
newDelta = np.vstack((Delta, h))
n, m = newL.shape[1], newL.shape[0]
transL = np.zeros((m*3, n*3))
newL = newL.tolist()
for i in range(m):
    for j in range(n):
        transL[i*3][j*3] = transL[i*3+1][j*3+1] = transL[i*3+2][j*3+2] = newL[i][j]
transL = scipy.sparse.coo_matrix(transL)
newDelta = newDelta.tolist()
transDelta = []
for item in newDelta:
    for x in item:
        transDelta.append(x)
transDelta=np.array(transDelta)
X = scipy.sparse.linalg.lsqr(transL, transDelta)[0]

# display
newmesh = mesh.copy()
for i in range(n):
    newmesh.vertices[old_id[i]] = [X[i*3],X[i*3+1],X[i*3+2]]
scene2 = [newmesh]
trimesh.Scene(scene2).show(viewer='gl',smooth=False)

    

成果：