Make a rest_post yourself. It is an equilateral triangle. The poses are a right triangle and an obtuse triangle.
You can run it directly to observe the effect
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import lsq_linear
from scipy.cluster.vq import vq, kmeans, whiten
import time
#PointCloudRegistrationAlgorithm
import numpy as np
def kabsch_2d(P, Q):
"""
Two-dimensional point cloud registration algorithm, using Kabsch algorithm.
inputs: P N x 2 numpy matrix, representing the two-dimensional coordinate points in the point set P
Q N x 2 numpy matrix, representing the two-dimensional coordinate points in the point set Q
return: A 3 x 2 matrix, the first two lines are rotation matrices, and the last line is translation vector
"""
if (P.size == 0 or Q.size == 0):
raise ValueError("Empty matrices sent to kabsch")
centroid_P = np.mean(P, axis=0)
centroid_Q = np.mean(Q, axis=0)
P_centered = P - centroid_P
Q_centered = Q - centroid_Q
H = P_centered.T.dot(Q_centered)
U, S, V = np.linalg.svd(H)
R = U.dot(V).T
if np.linalg.det(R) < 0:
V[1, :] *= -1
R = U.dot(V).T
t = centroid_Q - R.dot(centroid_P)
return np.vstack((R, t))
def initialize_2d(poses, rest_pose, num_bones, iterations=5):
"""
Uses the k-means algorithm to initialize bone transformations.
inputs: poses |num_poses| x |num_verts| x 2 matrix representing coordinates of vertices of each pose
rest_pose |num_verts| x 2 numpy matrix representing the coordinates of vertices in rest pose
num_bones Number of bones to initialize
Number of iterations to run the k-means algorithm
return: A |num_bones| x |num_poses| x 3 x 2 matrix representing the stacked Rotation and Translation
for each pose, for each bone.
A |num_bones| x 2 matrix representing the translations of the rest bones.
"""
#rest_pose: Modeling pose
num_verts = rest_pose.shape[0]
num_poses = poses.shape[0]
# trans for each bone per frame
# bone_transforms: (num_bones, num_poses, 3, 2)
bone_transforms = np.empty((num_bones, num_poses, 3, 2)) # [(R, T) for for each pose] for each bone
# 3rd dim has 3 rows for R and 1 row for T
rest_bones_t = np.empty((num_bones, 2)) # Translations for bones at rest pose
rest_pose_corrected = np.empty((num_bones, num_verts, 2)) # Rest pose - mean of vertices attached to each bone
# Use k-means to assign bones to vertices
# Whitening: rest_pose / The standard deviation of each column, that is, the standard deviation of the three axes of x, y, and z. This is only used for kmeans
whitened = whiten(rest_pose)
# Directly group the whitened rest poses into kmeans groups. The number of groups is the number of bones.
# python uses two functions to complete this step
# kmeans returns the center point
codebook, _ = kmeans(whitened, num_bones)#Clustering center (initial bone point)
# Bind the whitened initial vertex to the center point to complete the grouping
#vert_assignments: (num_verts,)
vert_assignments, _ = vq(whitened, codebook) #The index of the bone point corresponding to each vertex [0,1,5,6,8,.....]
# Traverse all bones
# Compute initial random bone transformations
for bone in range(num_bones):
# Calculate the center point. Each vert has a bone assignment. All vertices of each bone are used as a cluster to calculate the center.
# In rest pose, the current bone corresponds to the center of the cluster
#rest_bones_t:(num_bones, 2)
rest_bones_t[bone] = np.mean(rest_pose[vert_assignments == bone], axis=0)
# This corrected is the offset of each vertex relative to the center point of the cluster it belongs to.
# rest_pose: (num_verts, 3)
# rest_pose_corrected: (num_bones, num_verts, 3)
# Subtract the center point from each vertex in the rest pose. What needs to be noted here is that here are all the calculated vertices.
# In other words, here is the position of all the vertices of the body relative to the center of the cluster of a bone.
rest_pose_corrected[bone] = rest_pose - rest_bones_t[bone]
#rest_pose_corrected[bone] = rest_pose - np.mean(rest_pose[vert_assignments == bone], axis=0)
# Now we have the rest pose and then iterate through each pose
for pose in range(num_poses):
# The coordinates of the rest pose passed in relative to the center point, and the vertex data of the pose in the current frame
# For a bone, there may be multiple vertices bound to this bone, so for a bone, it corresponds to a point set
# Therefore, the goal of the kabsch algorithm is to try to rotate and translate this bone to minimize the difference between the two point sets.
# Calculate the transform of the bone in the current frame
# rest_pose_corrected contains all the vertices, but here only the vertices belonging to this bone are taken, and the kabsch algorithm is used to calculate the transform of the current bone.
bone_transforms[bone, pose] = kabsch_2d(rest_pose_corrected[bone, vert_assignments == bone], poses[pose, vert_assignments == bone])
for it in range(iterations):
# Re-assign bones to vertices using smallest reconstruction error from all poses
# This reconstructs the space of each bone, each frame, and all vertices.
# constructed: (num_bones, num_poses, num_verts, 3)
constructed = np.empty((num_bones, num_poses, num_verts, 2)) # |num_bones| x |num_poses| x |num_verts| x 3
for bone in range(num_bones):
# Traverse each bone
#Transform vertices into model space
# Rotate first
# Since the matrix has 4 rows and 3 columns, the first three rows of the matrix are taken here, which is the rotation matrix.
# Set the vertex position of the rest pose to the corresponding bone space, and then rotate it
# After the vertex is rotated, the yz axes of the vertex are reversed and become right-handed. The kabsch function constructing bone_transforms above also has the operation of transforming the hand system.
# But this is equivalent to all vertices in rest pose being transformed once in all frames.
# bone_transforms: (num_bones, num_poses, 4, 3)
Rp = bone_transforms[bone,:,:2,:].dot((rest_pose - rest_bones_t[bone]).T).transpose((0, 2, 1)) # |num_poses| x |num_verts| x 3
#Move again
# R * p + T
constructed[bone] = Rp + bone_transforms[bone, :, np.newaxis, 2, :]
# Compare the transformed vertex position with the original pose
# Calculate the 2 norm, that is, the sum of squares and then the square root
# The following axis=(1,3), 1 corresponds to poses, 3 corresponds to position, that is, the 2 norm is calculated based on the difference in vertex positions of all frames.
# That is, traverse num_bones and num_verts, sum the squares of all elements and then take the root sign
# constructed: (num_bones, num_poses, num_verts, 3)
# errs: (num_bones, num_verts)
errs = np.linalg.norm(constructed - poses, axis=(1, 3))
# print(errs.sum())
# Above we calculated errors. Its structure is (num_bones, num_verts), which is the error of each bone and each vertex (regardless of whether it belongs to the bone) in the entire animation sequence.
# For each bone, search for the smallest vertex index by column, that is, the index of the bone with the smallest error is stored at the position of each vertex. So we have a new mapping
# (num_verts),
vert_assignments = np.argmin(errs, axis=0)
## Visualization of vertex assignments for bone 0 over iterations
## Make 5 copies of an example pose mesh and call them test0, test1...
#for i in range(num_verts):
# if vert_assignments[i] == 0:
# pm.select('test{0}.vtx[{1}]'.format(it, i), add=True)
#print(vert_assignments)
# For each bone, for each pose, compute new transform using kabsch
for bone in range(num_bones):
# Recalculate the cluster center point based on the new bone mapping
rest_bones_t[bone] = np.mean(rest_pose[vert_assignments == bone], axis=0)
# Calculate the position of all vertices relative to each bone
rest_pose_corrected[bone] = rest_pose - np.mean(rest_pose[vert_assignments == bone], axis=0)
# Reuse kabsch algorithm to calculate bone transformation
for pose in range(num_poses):
# This is still based on vert_assignments and all vertices under a certain bone to calculate the transform of the bone.
bone_transforms[bone, pose] = kabsch_2d(rest_pose_corrected[bone, vert_assignments == bone], poses[pose, vert_assignments == bone])
#The final return is the animation data of the skeleton and the position of the skeleton in rest pose
return bone_transforms, rest_bones_t
def update_weight_map_2d(bone_transforms, rest_bones_t, poses, rest_pose, sparseness):
"""
Update the bone-vertex weight map W by fixing bone transformations and using a least squares
solver subject to non-negativity constraint, affinity constraint, and sparseness constraint.
inputs: bone_transforms |num_bones| x |num_poses| x 4 x 3 matrix representing the stacked
Rotation and Translation for each pose, for each bone.
rest_bones_t |num_bones| x 3 matrix representing the translations of the rest bones
poses |num_poses| x |num_verts| x 3 matrix representing coordinates of vertices of each pose
rest_pose |num_verts| x 3 numpy matrix representing the coordinates of vertices in rest pose
sparseness Maximum number of bones allowed to influence a particular vertex
return: A |num_verts| x |num_bones| weight map representing the influence of the jth bone on the ith vertex
"""
num_verts = rest_pose.shape[0]
num_poses = poses.shape[0]
num_bones = bone_transforms.shape[0]
W = np.empty((num_verts, num_bones))
for v in range(num_verts):
# Traverse all vertices
# For every vertex, solve a least squares problem
Rp = np.empty((num_bones, num_poses, 2))
# Transform the current vertex, use all bones to transform, and get the transformed vertex
for bone in range(num_bones):
Rp[bone] = bone_transforms[bone,:,:2,:].dot(rest_pose[v] - rest_bones_t[bone]) # |num_bones| x |num_poses| x 3
# R * p + T
Rp_T = Rp + bone_transforms[:, :, 2, :] # |num_bones| x |num_poses| x 3
# This is the transformed vertex position. This is equivalent to expanding each xyz of all poses by rows.
# (num_bones, nun_poses, 3)->(num_poses, 3, num_bones) -> (3 * num_poses, num_bones)
A = Rp_T.transpose((1, 2, 0)).reshape((2 * num_poses, num_bones)) # 3 * |num_poses| x |num_bones|
# This is the real vertex position
b = poses[:, v, :].reshape(2 * num_poses) # 3 * |num_poses| x 1
# Calculate when x has the smallest difference between Ax - b, using mse. Linear regression with bound is used here, because the values of all weights can only be between [0], 1]
#Finally got the weight of each bone
# Bounds ensure non-negativity constraint and kind of affinity constraint
w = lsq_linear(A, b, bounds=(0, 1), method='bvls').x # |num_bones| x 1
# Bone weight normalization
w /= np.sum(w) # Ensure that w sums to 1 (affinity constraint)
# Since each vertex has a maximum weight limit |K|, we have to discard other
# Remove |B| - |K| bone weights with the least "effect"
# Reshape A*w here and get the fitted vertices.
# (3 * num_poses, num_bones) -> (num_poses, 3, num_bones)
# Then based on column 1, find the second norm, that is, traverse num_poses and num_bones, and sum up the square root of all elements, which is actually the 2 norm of the vertex position.
effect = np.linalg.norm((A * w).reshape(num_poses, 2, num_bones), axis=1) # |num_poses| x |num_bones|
# Then merged num_poses. What we actually do here is to find out the degree to which the current vertex is affected by different bones. Then we can know which bones should be killed.
effect = np.sum(effect, axis=0) # |num_bones| x 1
# Because a vertex can only be affected by |k| root bones at most, the number of bones to be discarded is calculated here.
num_discarded = max(num_bones - sparseness, 0)
# np.argpartition is to group the contents of the array into groups. Those smaller than num_discarded are placed in front, and those larger than num_discarded are placed in the back.
# So here is actually discarding the smallest num_discarded elements
# The effective here saves the index of the selected bone for each vertex: the number is: |sparseness|
effective = np.argpartition(effect, num_discarded)[num_discarded:] # |sparseness| x 1
# At this point, we know which bones are bound to each vertex.
# Only use the few bones we selected and re-linear regression to calculate the weights
# Run least squares again, but only use the most effective bones
A_reduced = A[:, effective] # 3 * |num_poses| x |sparseness|
w_reduced = lsq_linear(A_reduced, b, bounds=(0, 1), method='bvls').x # |sparseness| x 1
w_reduced /= np.sum(w_reduced) # Ensure that w sums to 1 (affinity constraint)
# Finally, the weight sparse matrix of the current vertex is obtained
w_sparse = np.zeros(num_bones)
# The bones we selected are set to the recalculated weights
w_sparse[effective] = w_reduced
w_sparse /= np.sum(w_sparse) # Ensure that w_sparse sums to 1 (affinity constraint)
W[v] = w_sparse
# Return the weight assignment relationship of all vertices
# (num_verts, num_bones)
return W
def update_bone_transforms_2d(W, bone_transforms, rest_bones_t, poses, rest_pose):
"""
Updates the bone transformations by fixing the bone-vertex weight map and minimizing an
objective function individually for each pose and each bone.
inputs: W |num_verts| x |num_bones| matrix: bone-vertex weight map. Rows sum to 1, sparse.
bone_transforms |num_bones| x |num_poses| x 4 x 3 matrix representing the stacked
Rotation and Translation for each pose, for each bone.
rest_bones_t |num_bones| x 3 matrix representing the translations of the rest bones
poses |num_poses| x |num_verts| x 3 matrix representing coordinates of vertices of each pose
rest_pose |num_verts| x 3 numpy matrix representing the coordinates of vertices in rest pose
return: |num_bones| x |num_poses| x 4 x 3 matrix representing the stacked
Rotation and Translation for each pose, for each bone.
"""
num_bones = W.shape[1]
num_poses = poses.shape[0]
num_verts = W.shape[0]
# Traverse all poses
for pose in range(num_poses):
# Traverse all bones under the current pose
for bone in range(num_bones):
# Calculate the offset of all vertices relative to the current bone in rest pose
# Represents the points in rest pose without this rest bone's translation
p_corrected = rest_pose - rest_bones_t[bone] # |num_verts| x 3
# Calculate q_i for all vertices by equation (6)
# (num_bones, num_verts, 3)
constructed = np.empty((num_bones, num_verts, 2)) # |num_bones| x |num_verts| x 3
for bone2 in range(num_bones):
# Traverse all bones
# First calculate the rest pose based on bone2
# Extract the rotation matrix of the current bone, calculate the local position relative to the rest pose relative to the current bone, and rotate it
# can't use p_corrected before because we want to correct for every bone2 distinctly
# (num_verts, 3)
Rp = bone_transforms[bone2,pose,:2,:].dot((rest_pose - rest_bones_t[bone2]).T).T # |num_verts| x 3
# Then translate, here you get the vertex position transformed by each bone for each frame of animation.
# Since there are weights to indicate which vertices will be affected by which bones, for any bone, all vertices are counted here.
# R * p + T
constructed[bone2] = Rp + bone_transforms[bone2, pose, 2, :]
# w * (R * p + T)
# Here, multiply the vertices by the corresponding weights. Note that after the previous step, most of the weights here are 0. Only important weights have numerical values, and the maximum number is |K|
# So what is obtained here is the weighted position of each vertex relative to its corresponding bone. The sum is then required to obtain the true position.
# (num_bones, num_verts, 3)->(num_verts, num_bones, 3)
constructed = constructed.transpose((1, 0, 2)) * W[:, :, np.newaxis] # |num_verts| x |num_bones| x 3
# The operation here is to remove these positions from the current bone column to remove the influence of the current bone.
constructed = np.delete(constructed, bone, axis=1) # |num_verts| x |num_bones-1| x 3
# The current bone is the bone we are paying attention to, so q is the real vertex position minus the influence of other bones on the vertex, and the resulting point cloud
# Intuitive understanding is that pose is a real point cloud, and the constructed points without the current skeleton are an estimated point cloud.
q = poses[pose] - np.sum(constructed, axis=1) # |num_verts| x 3
# And p is the point cloud of the offset of the vertex in the rest pose relative to the center point. This center point is a weighted average, and the weight is related to the w of the vertex to the bone.
# p_star is the center point under rest pose
# q_star is the center point of the residual of the current animation
# Calculate p_star, q_star, p_bar, and q_bar for all verts by equation (8)
p_star = np.sum(np.square(W[:, bone, np.newaxis]) * p_corrected, axis=0) # |num_verts| x 3 => 3 x 1
p_star /= np.sum(np.square(W[:, bone])) # 3 x 1
q_star = np.sum(W[:, bone, np.newaxis] * q, axis=0) # |num_verts| x 3 => 3 x 1
q_star /= np.sum(np.square(W[:, bone])) # 3 x 1
# So far we have obtained two point clouds,
p_bar = p_corrected - p_star # |num_verts| x 3
q_bar = q - W[:, bone, np.newaxis] * q_star # |num_verts| x 3
# Next, perform point cloud registration to obtain the rotation and translation of the bones.
# Perform SVD by equation (9)
P = (p_bar * W[:, bone, np.newaxis]).T # 3 x |num_verts|
Q = q_bar.T # 3 x |num_verts|
# Point cloud registration to obtain rotation matrix
U, S, V = np.linalg.svd(np.matmul(P, Q.T))
# Calculate rotation R and translation t by equation (10)
R = U.dot(V).T # 3 x 3
t = q_star - R.dot(p_star) # 3 x 1
bone_transforms[bone, pose, :2, :] = R
bone_transforms[bone, pose, 2, :] = t
# Finally return, the pose of each bone in each frame
return bone_transforms
def reconstruction_err_2d(poses, rest_pose, bone_transforms, rest_bones_t, W):
num_bones = bone_transforms.shape[0]
num_verts = W.shape[0]
num_poses = poses.shape[0]
# Calculate the local offset of the rest pose relative to each bone
p_corrected = rest_pose[np.newaxis, :, :] - rest_bones_t[:, np.newaxis, :] # |num_bones| x |num_verts| x 2
constructions = np.empty((num_bones, num_poses, num_verts, 2)) # |num_bones| x |num_poses| x |num_verts| x 2
#Perform rotation transformation
for bone in range(num_bones):
constructions[bone] = np.einsum('ijk,lk->ilj', bone_transforms[bone, :, :2, :], p_corrected[bone]) # |num_poses| x |num_verts| x 2
#Add offset
constructions + = bone_transforms[:, :, np.newaxis, 2, :] # |num_bones| x |num_poses| x |num_verts| x 2
# Calculate the weights to get the final fitting vertex
constructions *= (W.T)[:, np.newaxis, :, np.newaxis] # |num_bones| x |num_poses| x |num_verts| x 2
# Subtract the fitted pose from the real pose to get the residual errors
constructions = np.sum(constructions, axis=0)
errors = poses - constructions # |num_poses| x |num_verts| x 2
#draw图
plt.figure(figsize=(6, 6))
plt.scatter(constructions[1, :, 0], constructions[1, :, 1], s=5, color='blue')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Points Inside Equilateral Triangle')
plt.xlim(1, 2)
plt.ylim(1, 2)
plt.grid(True)
plt.show()
# Calculate the second norm as the final error
return np.mean(np.linalg.norm(errors, axis=2))
def SSDR_2d(poses, rest_pose, num_bones, sparseness=4, max_iterations=20):
"""
Computes the Smooth Skinning Decomposition with Rigid bones
inputs: poses |num_poses| x |num_verts| x 2 matrix representing coordinates of vertices of each pose
rest_pose |num_verts| x 2 numpy matrix representing the coordinates of vertices in rest pose
num_bones number of bones to create
sparseness max number of bones influencing a single vertex
return: An i x j matrix of bone-vertex weights, where i = # vertices and j = # bones
A length-B list of (length-t lists of bone transformations [R_j | T_j] ), one list for each bone
A list of bone translations for the bones at rest
"""
start_time = time.time()
# Initialization, given animation frame poses, t pose, and how many bones are needed, return the transformation matrix of all bones and the center point of each bone
bone_transforms, rest_bones_t = initialize_2d(poses, rest_pose, num_bones)
print("init success")
for _ in range(max_iterations):
W = update_weight_map_2d(bone_transforms, rest_bones_t, poses, rest_pose, sparseness)
bone_transforms = update_bone_transforms_2d(W, bone_transforms, rest_bones_t, poses, rest_pose)
print("Reconstruction error:", reconstruction_err_2d(poses, rest_pose, bone_transforms, rest_bones_t, W))
end_time = time.time()
print("Done. Calculation took {0} seconds".format(end_time - start_time))
print("Avg reconstruction error:", reconstruction_err_2d(poses, rest_pose, bone_transforms, rest_bones_t, W))
return W, bone_transforms, rest_bones_t
if __name__ == '__main__':
# Define the vertices of an equilateral triangle
vertex1_equilateral = np.array([0, 0])
vertex2_equilateral = np.array([1, 0])
vertex3_equilateral = np.array([0.5, np.sqrt(3) / 2])
# Define the vertices of the right triangle
vertex1_right = np.array([0, 0])
vertex2_right = np.array([0, 1])
vertex3_right = np.array([1, 0])
# Define the vertices of an obtuse triangle
vertex1_blunt = np.array([0, 0])
vertex2_blunt = np.array([1, 0])
vertex3_blunt = np.array([0.2, 0.8])
# Generate 500 equally spaced points inside the triangle (equilateral triangle)
num_points = 40
t = np.linspace(0, 1, num_points)
points_inside_equilateral = np.array(
[(1 - a - b) * vertex1_equilateral + a * vertex2_equilateral + b * vertex3_equilateral
for a in t for b in t if a + b <= 1])
points_inside_equilateral = points_inside_equilateral.astype(float)
# Generate 500 equally spaced points inside the triangle (right triangle)
points_inside_right_triangle = np.array([(1 - a - b) * vertex1_right + a * vertex2_right + b * vertex3_right
for a in t for b in t if a + b <= 1])
points_inside_right_triangl = points_inside_right_triangle.astype(float)
# Generate 500 equally spaced points inside the triangle (obtuse triangle)
points_inside_blunt_triangle = np.array([(1 - a - b) * vertex1_blunt + a * vertex2_blunt + b * vertex3_blunt
for a in t for b in t if a + b <= 1])
points_inside_blunt_triangle= points_inside_blunt_triangle.astype(float)
#dynamic action
poses = np.zeros((2, 820, 2))
poses[0] = points_inside_right_triangle + 1
poses[1]=points_inside_blunt_triangle + 1
#static action
rest_pose = points_inside_equilateral + 1
SSDR_2d(poses,rest_pose,10,8,25)
This image is driven by equilateral triangles
The knowledge points of the article match the official knowledge files, and you can further learn related knowledge. Algorithm skill tree Home page Overview 57032 people are learning the system