iti/HandMesh
0
1import math2import heapq3import numpy as np4import os5import scipy.sparse as sp6from psbody.mesh import Mesh7 8 9def row(A):10 return A.reshape((1, -1))11 12 13def col(A):14 return A.reshape((-1, 1))15 16 17def get_vert_connectivity(mesh_v, mesh_f):18 """Returns a sparse matrix (of size #verts x #verts) where each nonzero19 element indicates a neighborhood relation. For example, if there is a20 nonzero element in position (15,12), that means vertex 15 is connected21 by an edge to vertex 12."""22 23 vpv = sp.csc_matrix((len(mesh_v), len(mesh_v)))24 25 # for each column in the faces...26 for i in range(3):27 IS = mesh_f[:, i]28 JS = mesh_f[:, (i + 1) % 3]29 data = np.ones(len(IS))30 ij = np.vstack((row(IS.ravel()), row(JS.ravel())))31 mtx = sp.csc_matrix((data, ij), shape=vpv.shape)32 vpv = vpv + mtx + mtx.T33 34 return vpv35 36 37def get_vertices_per_edge(mesh_v, mesh_f):38 """Returns an Ex2 array of adjacencies between vertices, where39 each element in the array is a vertex index. Each edge is included40 only once. If output of get_faces_per_edge is provided, this is used to41 avoid call to get_vert_connectivity()"""42 43 vc = sp.coo_matrix(get_vert_connectivity(mesh_v, mesh_f))44 result = np.hstack((col(vc.row), col(vc.col)))45 result = result[result[:, 0] < result[:, 1]] # for uniqueness46 47 return result48 49 50def vertex_quadrics(mesh):51 """Computes a quadric for each vertex in the Mesh.52 53 Returns:54 v_quadrics: an (N x 4 x 4) array, where N is # vertices.55 """56 57 # Allocate quadrics58 v_quadrics = np.zeros((59 len(mesh.v),60 4,61 4,62 ))63 64 # For each face...65 for f_idx in range(len(mesh.f)):66 67 # Compute normalized plane equation for that face68 vert_idxs = mesh.f[f_idx]69 verts = np.hstack((mesh.v[vert_idxs], np.array([1, 1,70 1]).reshape(-1, 1)))71 u, s, v = np.linalg.svd(verts)72 eq = v[-1, :].reshape(-1, 1)73 eq = eq / (np.linalg.norm(eq[0:3]))74 75 # Add the outer product of the plane equation to the76 # quadrics of the vertices for this face77 for k in range(3):78 v_quadrics[mesh.f[f_idx, k], :, :] += np.outer(eq, eq)79 80 return v_quadrics81 82 83def setup_deformation_transfer(source, target, use_normals=False):84 rows = np.zeros(3 * target.v.shape[0])85 cols = np.zeros(3 * target.v.shape[0])86 coeffs_v = np.zeros(3 * target.v.shape[0])87 coeffs_n = np.zeros(3 * target.v.shape[0])88 89 nearest_faces, nearest_parts, nearest_vertices = source.compute_aabb_tree(90 ).nearest(target.v, True)91 nearest_faces = nearest_faces.ravel().astype(np.int64)92 nearest_parts = nearest_parts.ravel().astype(np.int64)93 nearest_vertices = nearest_vertices.ravel()94 95 for i in range(target.v.shape[0]):96 # Closest triangle index97 f_id = nearest_faces[i]98 # Closest triangle vertex ids99 nearest_f = source.f[f_id]100 101 # Closest surface point102 nearest_v = nearest_vertices[3 * i:3 * i + 3]103 # Distance vector to the closest surface point104 dist_vec = target.v[i] - nearest_v105 106 rows[3 * i:3 * i + 3] = i * np.ones(3)107 cols[3 * i:3 * i + 3] = nearest_f108 109 n_id = nearest_parts[i]110 if n_id == 0:111 # Closest surface point in triangle112 A = np.vstack((source.v[nearest_f])).T113 coeffs_v[3 * i:3 * i + 3] = np.linalg.lstsq(A, nearest_v,114 rcond=-1)[0]115 elif n_id > 0 and n_id <= 3:116 # Closest surface point on edge117 A = np.vstack((source.v[nearest_f[n_id - 1]],118 source.v[nearest_f[n_id % 3]])).T119 tmp_coeffs = np.linalg.lstsq(A, target.v[i], rcond=-1)[0]120 coeffs_v[3 * i + n_id - 1] = tmp_coeffs[0]121 coeffs_v[3 * i + n_id % 3] = tmp_coeffs[1]122 else:123 # Closest surface point a vertex124 coeffs_v[3 * i + n_id - 4] = 1.0125 126 matrix = sp.csc_matrix((coeffs_v, (rows, cols)),127 shape=(target.v.shape[0], source.v.shape[0]))128 return matrix129 130 131def qslim_decimator_transformer(mesh, factor=None, n_verts_desired=None):132 """Return a simplified version of this mesh.133 134 A Qslim-style approach is used here.135 136 :param factor: fraction of the original vertices to retain137 :param n_verts_desired: number of the original vertices to retain138 :returns: new_faces: An Fx3 array of faces, mtx: Transformation matrix139 """140 141 if factor is None and n_verts_desired is None:142 raise Exception('Need either factor or n_verts_desired.')143 144 if n_verts_desired is None:145 n_verts_desired = math.ceil(len(mesh.v) * factor)146 147 Qv = vertex_quadrics(mesh)148 149 # fill out a sparse matrix indicating vertex-vertex adjacency150 # from psbody.mesh.topology.connectivity import get_vertices_per_edge151 vert_adj = get_vertices_per_edge(mesh.v, mesh.f)152 # vert_adj = sp.lil_matrix((len(mesh.v), len(mesh.v)))153 # for f_idx in range(len(mesh.f)):154 # vert_adj[mesh.f[f_idx], mesh.f[f_idx]] = 1155 156 vert_adj = sp.csc_matrix(157 (vert_adj[:, 0] * 0 + 1, (vert_adj[:, 0], vert_adj[:, 1])),158 shape=(len(mesh.v), len(mesh.v)))159 vert_adj = vert_adj + vert_adj.T160 vert_adj = vert_adj.tocoo()161 162 def collapse_cost(Qv, r, c, v):163 Qsum = Qv[r, :, :] + Qv[c, :, :]164 p1 = np.vstack((v[r].reshape(-1, 1), np.array([1]).reshape(-1, 1)))165 p2 = np.vstack((v[c].reshape(-1, 1), np.array([1]).reshape(-1, 1)))166 167 destroy_c_cost = p1.T.dot(Qsum).dot(p1)168 destroy_r_cost = p2.T.dot(Qsum).dot(p2)169 result = {170 'destroy_c_cost': destroy_c_cost,171 'destroy_r_cost': destroy_r_cost,172 'collapse_cost': min([destroy_c_cost, destroy_r_cost]),173 'Qsum': Qsum174 }175 return result176 177 # construct a queue of edges with costs178 queue = []179 for k in range(vert_adj.nnz):180 r = vert_adj.row[k]181 c = vert_adj.col[k]182 183 if r > c:184 continue185 186 cost = collapse_cost(Qv, r, c, mesh.v)['collapse_cost']187 heapq.heappush(queue, (cost, (r, c)))188 189 # decimate190 collapse_list = []191 nverts_total = len(mesh.v)192 faces = mesh.f.copy()193 while nverts_total > n_verts_desired:194 e = heapq.heappop(queue)195 r = e[1][0]196 c = e[1][1]197 if r == c:198 continue199 200 cost = collapse_cost(Qv, r, c, mesh.v)201 if cost['collapse_cost'] > e[0]:202 heapq.heappush(queue, (cost['collapse_cost'], e[1]))203 # print 'found outdated cost, %.2f < %.2f' % (e[0], cost['collapse_cost'])204 continue205 else:206 207 # update old vert idxs to new one,208 # in queue and in face list209 if cost['destroy_c_cost'] < cost['destroy_r_cost']:210 to_destroy = c211 to_keep = r212 else:213 to_destroy = r214 to_keep = c215 216 collapse_list.append([to_keep, to_destroy])217 218 # in our face array, replace "to_destroy" vertidx with "to_keep" vertidx219 np.place(faces, faces == to_destroy, to_keep)220 221 # same for queue222 which1 = [223 idx for idx in range(len(queue))224 if queue[idx][1][0] == to_destroy225 ]226 which2 = [227 idx for idx in range(len(queue))228 if queue[idx][1][1] == to_destroy229 ]230 for k in which1:231 queue[k] = (queue[k][0], (to_keep, queue[k][1][1]))232 for k in which2:233 queue[k] = (queue[k][0], (queue[k][1][0], to_keep))234 235 Qv[r, :, :] = cost['Qsum']236 Qv[c, :, :] = cost['Qsum']237 238 a = faces[:, 0] == faces[:, 1]239 b = faces[:, 1] == faces[:, 2]240 c = faces[:, 2] == faces[:, 0]241 242 # remove degenerate faces243 def logical_or3(x, y, z):244 return np.logical_or(x, np.logical_or(y, z))245 246 faces_to_keep = np.logical_not(logical_or3(a, b, c))247 faces = faces[faces_to_keep, :].copy()248 249 nverts_total = (len(np.unique(faces.flatten())))250 251 new_faces, mtx = _get_sparse_transform(faces, len(mesh.v))252 return new_faces, mtx253 254 255def _get_sparse_transform(faces, num_original_verts):256 verts_left = np.unique(faces.flatten())257 IS = np.arange(len(verts_left))258 JS = verts_left259 data = np.ones(len(JS))260 261 mp = np.arange(0, np.max(faces.flatten()) + 1)262 mp[JS] = IS263 new_faces = mp[faces.copy().flatten()].reshape((-1, 3))264 265 ij = np.vstack((IS.flatten(), JS.flatten()))266 mtx = sp.csc_matrix((data, ij),267 shape=(len(verts_left), num_original_verts))268 269 return (new_faces, mtx)270 271 272def generate_transform_matrices(mesh, factors):273 """Generates len(factors) meshes, each of them is scaled by factors[i] and274 computes the transformations between them.275 276 Returns:277 M: a set of meshes downsampled from mesh by a factor specified in factors.278 A: Adjacency matrix for each of the meshes279 D: csc_matrix Downsampling transforms between each of the meshes280 U: Upsampling transforms between each of the meshes281 F: a list of faces282 """283 284 factors = map(lambda x: 1.0 / x, factors)285 M, A, D, U, F, V = [], [], [], [], [], []286 F.append(mesh.f) # F[0]287 V.append(mesh.v)288 A.append(get_vert_connectivity(mesh.v, mesh.f).astype('float32')) # A[0]289 M.append(mesh) # M[0]290 291 for factor in factors:292 ds_f, ds_D = qslim_decimator_transformer(M[-1], factor=factor)293 D.append(ds_D.astype('float32'))294 new_mesh_v = ds_D.dot(M[-1].v)295 new_mesh = Mesh(v=new_mesh_v, f=ds_f)296 F.append(new_mesh.f)297 V.append(new_mesh.v)298 M.append(new_mesh)299 A.append(300 get_vert_connectivity(new_mesh.v, new_mesh.f).astype('float32'))301 U.append(setup_deformation_transfer(M[-1], M[-2]).astype('float32'))302 303 return M, A, D, U, F, V304