mesh subdivision (#13)

README.md

@@ -145,8 +145,7 @@ properly credited both in this page and in the individual files.
 - Vectorize `write_ply.py`
 - Fix number type thing in c++ bindings
 - Write `read_ply.py`
-- Write `subdivide_triangles.py` and `subdivide_tets.py`
-- Write `loop.py` for loop subdivision (returning sparse mapping matrix)
+- Add tets to `subdivide.py`
 - Write quadratic solver with fixed points and linear constraints
 - Write `decimate.py` functionality
 - Explore exactly which part of png2poly's dependencies we need

gpytoolbox/__init__.py

@@ -79,3 +79,4 @@
 from .triangle_triangle_adjacency import triangle_triangle_adjacency
 from .halfedge_edge_map import halfedge_edge_map
 from .array_correspondence import array_correspondence
+from .subdivide import subdivide

gpytoolbox/subdivide.py

@@ -0,0 +1,265 @@
+import numpy as np
+import scipy as sp
+from gpytoolbox.halfedge_edge_map import halfedge_edge_map
+
+def subdivide(V,F,
+    method='upsample',
+    iters=1,
+    return_matrix=False):
+    # Subdivides a mesh
+    #
+    # Input:
+    #       V  #V by d numpy array of mesh vertex positions.
+    #          Can be None, in which case the subdivision is performed only
+    #          logically on the faces.
+    #       F  #F by k int numpy array of vertex indices into V.
+    #          Can represent a polyline or a triangle mesh.
+    #       Optional:
+    #                method  which method to use for subdivison.
+    #                        can be 'upsample' {default},
+    #                               'loop' (only triangle meshes)
+    #                iters  how many subdivision iterations to perform
+    #                return_matrix  whether to return matrix for sparse map S
+    #
+    # Output:
+    #       Vu  #Vu by d numpy array of subdivided mesh vertex positions
+    #           Is None if V was None.
+    #       Fu  #Fu by k int numpy array of subdivided vertex indices into V.
+    #       S  If return_matrix is True, the sparse matrix such that Vu == S*V
+    #
+
+    assert iters >= 0
+
+    Vu,Fu = V,F
+    if return_matrix:
+        S = sp.sparse.eye(_n(V,F), format='csr')
+
+    for i in range(iters):
+        Vu,Fu,St = _one_subdivision(Vu,Fu,method,return_matrix)
+        if return_matrix:
+            S = St*S
+
+    # If there were no iterations at all, we want to return a copy of the
+    # original.
+    if iters==0:
+        Vu = None if V is None else V.copy()
+        Fu = F.copy()
+
+    if return_matrix:
+        return Vu, Fu, S
+    else:
+        return Vu, Fu
+
+
+def _one_subdivision(V,F,method,return_matrix):
+    # Dispatcher function for one single subdivision.
+    k = F.shape[1]
+    if k==2:
+        if method=='upsample':
+            Vu,Fu,S = _upsample_polyline(V,F,return_matrix)
+        else:
+            assert False, "Method not supported for this simplex size."
+    elif k==3:
+        if method=='upsample':
+            Vu,Fu,S = _upsample_triangle_mesh(V,F,return_matrix)
+        elif method=='loop':
+            Vu,Fu,S = _loop_triangle_mesh(V,F,return_matrix)
+        else:
+            assert False, "Method not supported for this simplex size."
+
+    else:
+        assert False, "Simplex dimension not supported."
+
+    return Vu, Fu, S
+
+
+def _upsample_polyline(V,F,return_matrix):
+    # Upsampling of every polyline by inserting a vertex in the middle of each
+    # segment.
+    assert V.shape[0]>1
+    assert F.shape[1]==2
+    assert F.shape[0]>0
+
+    n = _n(V,F)
+    m = F.shape[0]
+    new = np.arange(n,n+m)
+    Fu = np.block([
+        [F[:,0], new],
+        [new,    F[:,1]]
+        ]).transpose()
+
+    if V is None:
+        Vu = None
+    else:
+        Vnew = 0.5*(V[F[:,0],:] + V[F[:,1],])
+        Vu = np.concatenate((V,Vnew), axis=0)
+
+    if return_matrix:
+        old = np.arange(0,n)
+        i = np.concatenate((
+            old,
+            new,
+            new
+            ))
+        j = np.concatenate((
+            old,
+            F[:,0],
+            F[:,1]
+            ))
+        k = np.concatenate((
+            np.full(n,1.),
+            np.full(m,0.5),
+            np.full(m,0.5)
+            ))
+        S = sp.sparse.csr_matrix((k,(i,j)), shape=(n+m,n))
+    else:
+        S = None
+
+    return Vu,Fu,S
+
+
+def _upsample_triangle_mesh(V,F,return_matrix,
+    return_halfedge_edge_map=False):
+    # Upsampling of every triangle by inserting a vertex in the middle of each
+    # edge.
+    assert F.shape[1]==3
+    assert F.shape[0]>0
+
+    n = _n(V,F)
+    m = F.shape[0]
+    he,E,he_to_E,E_to_he = halfedge_edge_map(F, assume_manifold=True)
+    e = E.shape[0]
+
+    # A triangle [0,1,2] is replaced with four new triangles:
+    # [[0,3,5], [3,1,4], [5,4,2], [5,3,4]]
+    # Vertex k+3 is in the middle of the edge from k to k+1
+    Fu = np.block([
+        [F[:,0],         n+he_to_E[:,2], n+he_to_E[:,1], n+he_to_E[:,1]],
+        [n+he_to_E[:,2], F[:,1],         n+he_to_E[:,0], n+he_to_E[:,2]],
+        [n+he_to_E[:,1], n+he_to_E[:,0], F[:,2],         n+he_to_E[:,0]]
+        ]).transpose()
+
+    if V is None:
+        Vu = None
+    else:
+        assert V.shape[0]>1
+        Vnew = 0.5*(V[E[:,0],:] + V[E[:,1],:])
+        Vu = np.concatenate((V,Vnew), axis=0)
+
+    if return_matrix:
+        old = np.arange(0,n)
+        new = np.arange(n,n+e)
+        i = np.concatenate((
+            old,
+            new,
+            new
+            ))
+        j = np.concatenate((
+            old,
+            E[:,0],
+            E[:,1]
+            ))
+        k = np.concatenate((
+            np.full(n,1.),
+            np.full(e,0.5),
+            np.full(e,0.5)
+            ))
+        S = sp.sparse.csr_matrix((k,(i,j)), shape=(n+e,n))
+    else:
+        S = None
+
+    if return_halfedge_edge_map:
+        return Vu,Fu,S,he,E,he_to_E,E_to_he
+    else:
+        return Vu,Fu,S
+
+def _loop_triangle_mesh(V,F,return_matrix):
+    # Apply one iteration of Loop subdivision.
+    # This, by default, uses the movement rule and β from
+    # https://graphics.stanford.edu/~mdfisher/subdivision.html
+    # (not Loop's original β)
+
+    n = _n(V,F)
+    m = F.shape[0]
+    _,Fu,_,he,E,he_to_E,E_to_he = \
+    _upsample_triangle_mesh(None,F,False,return_halfedge_edge_map=True)
+    e = E.shape[0]
+
+    # Group halfedges
+    bE_mask = (E_to_he[:,1,:] == -1).any(axis=-1)
+    bdry_he = E_to_he[bE_mask,0,:]
+    int_he = np.concatenate((E_to_he[~bE_mask,0,:],E_to_he[~bE_mask,1,:]), axis=0)
+    iV_mask = np.full(n, True)
+    iV_mask[E[bE_mask,:].ravel()] = False
+    int_tail_he = np.stack(np.nonzero(iV_mask[he[:,:,0]]), axis=-1)
+
+    # Compute, for each vertex, the β needed for Loop.
+    n_adj = np.bincount(E.ravel(), minlength=n)
+    n_adj[n_adj==0] = -1
+    β = np.where(n_adj<3, np.NAN, np.where(n_adj==3, 3./16., (3./8.)/n_adj))
+
+    # We always compute the matrix S since we need it to construct Vu
+    i = np.concatenate((
+        # Add new boundary vertex in the middle of each bdry_he halfedge
+        n+he_to_E[bdry_he[:,0], bdry_he[:,1]],
+        n+he_to_E[bdry_he[:,0], bdry_he[:,1]],
+        # Move old bdry vertex based on 3/8 1/8 1/8 rule
+        he[bdry_he[:,0], bdry_he[:,1], 0],
+        he[bdry_he[:,0], bdry_he[:,1], 0],
+        he[bdry_he[:,0], bdry_he[:,1], 1],
+        # Add new interior vertex with 3/8 3/8 1/8 1/8 rule
+        n+he_to_E[int_he[:,0], int_he[:,1]],
+        n+he_to_E[int_he[:,0], int_he[:,1]],
+        # Move old interior vertex with β, 1-βn rule
+        he[int_tail_he[:,0], int_tail_he[:,1], 0],
+        he[int_tail_he[:,0], int_tail_he[:,1], 0]
+        ))
+    j = np.concatenate((
+        # Add new boundary vertex in the middle of each bdry_he halfedge
+        he[bdry_he[:,0], bdry_he[:,1], 0],
+        he[bdry_he[:,0], bdry_he[:,1], 1],
+        # Move old bdry vertex based on 3/8 1/8 1/8 rule
+        he[bdry_he[:,0], bdry_he[:,1], 0],
+        he[bdry_he[:,0], bdry_he[:,1], 1],
+        he[bdry_he[:,0], bdry_he[:,1], 0],
+        # Add new interior vertex with 3/8 3/8 1/8 1/8 rule
+        he[int_he[:,0], int_he[:,1], 0],
+        F[int_he[:,0], int_he[:,1]],
+        # Move old interior vertex with β, 1-βn rule
+        he[int_tail_he[:,0], int_tail_he[:,1], 0],
+        he[int_tail_he[:,0], int_tail_he[:,1], 1]
+        ))
+    k = np.concatenate((
+        # Add new boundary vertex in the middle of each bdry_he halfedge
+        np.full(bdry_he.shape[0], 0.5),
+        np.full(bdry_he.shape[0], 0.5),
+        # Move old bdry vertex based on 3/4 1/8 1/8 rule
+        np.full(bdry_he.shape[0], 0.75),
+        np.full(bdry_he.shape[0], 0.125),
+        np.full(bdry_he.shape[0], 0.125),
+        # Add new interior vertex with 3/8 3/8 1/8 1/8 rule
+        np.full(int_he.shape[0], 0.375),
+        np.full(int_he.shape[0], 0.125),
+        # Move old interior vertex with β, 1-βn rule
+        (1./n_adj[he[int_tail_he[:,0], int_tail_he[:,1],0]] - 
+            β[he[int_tail_he[:,0] , int_tail_he[:,1],0]]),
+        β[he[int_tail_he[:,0], int_tail_he[:,1],0]]
+        ))
+    S = sp.sparse.csr_matrix((k,(i,j)), shape=(n+e,n))
+
+    if V is None:
+        Vu = None
+    else:
+        Vu = S*V
+
+    if not return_matrix:
+        S = None
+
+    return Vu,Fu,S
+
+
+def _n(V,F):
+    if V is None:
+        return np.max(F) + 1
+    else:
+        return V.shape[0]

test/test_halfedge_edge_map.py

@@ -70,7 +70,18 @@ def variety_of_asserts(self, he, E, he_to_E, E_to_he, assume_manifold):
                 self.assertTrue(np.all(np.sort(he[he_ind_valid], axis=-1) == np.sort(E[e,:], axis=-1)))
         else:
             self.assertTrue(all([np.all(np.sort(he[E_to_he[e][:,0],E_to_he[e][:,1],:], axis=-1) == np.sort(E[e,:], axis=-1)) for e in range(len(E_to_he))]))
-
+        # Test reversibility
+        for i in range(he_to_E.shape[0]):
+            for j in range(3):
+                e = he_to_E[i,j]
+                if (E_to_he[e][0,:]==np.array([i,j])).all():
+                    self.assertTrue(E_to_he[e].shape==(1,2) or
+                        (E_to_he[e][1,:]==np.array([-1,-1])).all() or
+                        he_to_E[E_to_he[e][1,0],E_to_he[e][1,1]]==e)
+                else:
+                    self.assertTrue((E_to_he[e][1,:]==np.array([i,j])).all())
+                    self.assertTrue((E_to_he[e][0,:]==np.array([-1,-1])).all() or
+                        he_to_E[E_to_he[e][0,0],E_to_he[e][0,1]]==e)
 
 if __name__ == '__main__':
     unittest.main()