pygel3d.hmesh

  1""" The hmesh module provides an halfedge based mesh representation. 
  2In addition this module contains a variety of functions for mesh manipulation and inspection.
  3Specifcally, the module contains functions for mesh simplification, smoothing, subdivision, 
  4and editing of vertices, faces, and edges. The volumetric_isocontour function allows us to 
  5create a polygonal mesh from volumetric data by isocontouring. The skeleton_to_feq function 
  6allows us to turn a skeleton graph into a Face Extrusion Quad Mesh.
  7"""
  8import ctypes as ct
  9import numpy as np
 10from math import sqrt
 11from numpy import ndarray
 12from numpy.linalg import norm
 13from pygel3d import lib_py_gel, IntVector, Vec3dVector, spatial
 14from scipy.sparse import csc_matrix, vstack
 15from scipy.sparse.linalg import lsqr
 16from collections import defaultdict
 17
 18
 19class Manifold:
 20    """ The Manifold class represents a halfedge based mesh. It is maybe a bit grand to call
 21    a mesh class Manifold, but meshes based on the halfedge representation are manifold (if we
 22    ignore a few corner cases) unlike some other representations. This class contains a number of
 23    methods for mesh manipulation and inspection. Note also that numerous further functions are
 24    available to manipulate meshes stored as Manifolds.
 25
 26    Many of the functions below accept arguments called hid, fid, or vid. These are simply indices
 27    of halfedges, faces and vertices, respectively: integer numbers that identify the corresponding
 28    mesh element. Using a plain integer to identify a mesh entity means that, for instance, a
 29    vertex index can also be used as an index into, say, a NumPy array without any conversion.
 30    """
 31    def __init__(self,orig=None):
 32        if orig == None:
 33            self.obj = lib_py_gel.Manifold_new()
 34        else:
 35            self.obj = lib_py_gel.Manifold_copy(orig.obj)
 36    @classmethod
 37    def from_triangles(cls,vertices, faces):
 38        """ Given a list of vertices and triangles (faces), this function produces
 39        a Manifold mesh."""
 40        m = cls()
 41        vertices = np.asarray(vertices,dtype=np.float64, order='C')
 42        faces = np.asarray(faces,dtype=ct.c_int, order='C')
 43        m.obj = lib_py_gel.Manifold_from_triangles(len(vertices),len(faces),vertices,faces)
 44        return m
 45    @classmethod
 46    def from_points(cls,pts,xaxis=np.array([1,0,0]),yaxis=np.array([0,1,0])):
 47        """ This function computes the Delaunay triangulation of pts. You need
 48        to specify xaxis and yaxis if they are not canonical. The function returns
 49        a Manifold with the resulting triangles. Clearly, this function will
 50        give surprising results if the surface represented by the points is not
 51        well represented as a 2.5D surface, aka a height field. """
 52        m = cls()
 53        pts = np.asarray(pts,dtype=np.float64, order='C')
 54        xaxis = np.asarray(xaxis,dtype=np.float64, order='C')
 55        yaxis = np.asarray(yaxis,dtype=np.float64, order='C')
 56        m.obj = lib_py_gel.Manifold_from_points(len(pts), pts, xaxis, yaxis)
 57        return m
 58    def __del__(self):
 59        lib_py_gel.Manifold_delete(self.obj)
 60    def merge_with(self, other):
 61        """ Merge this Manifold with another one given as the argument. This function
 62        does not return anything. It simply combines the two meshes in the Manifold on which
 63        the method is called. """
 64        lib_py_gel.Manifold_merge(self.obj, other.obj)
 65    def add_face(self,pts):
 66        """ Add a face to the Manifold.
 67        This function takes a list of 3D points, pts, as argument and creates a face
 68        in the mesh with those points as vertices. The function returns the index
 69        of the created face.
 70        """
 71        pts = np.asarray(pts,dtype=np.float64, order='C')
 72        return lib_py_gel.Manifold_add_face(self.obj, len(pts), pts)
 73    def positions(self):
 74        """ Retrieve an array containing the vertex positions of the Manifold.
 75        It is not a copy: any changes are made to the actual vertex positions. """
 76        pos = ct.POINTER(ct.c_double)()
 77        n = lib_py_gel.Manifold_positions(self.obj, ct.byref(pos))
 78        return np.ctypeslib.as_array(pos,(n,3))
 79    def no_allocated_vertices(self):
 80        """ Number of vertices.
 81        This number could be higher than the number of actually
 82        used vertices, but corresponds to the size of the array allocated
 83        for vertices."""
 84        return lib_py_gel.Manifold_no_allocated_vertices(self.obj)
 85    def no_allocated_faces(self):
 86        """ Number of faces.
 87        This number could be higher than the number of actually
 88        used faces, but corresponds to the size of the array allocated
 89        for faces."""
 90        return lib_py_gel.Manifold_no_allocated_faces(self.obj)
 91    def no_allocated_halfedges(self):
 92        """ Number of halfedges.
 93        This number could be higher than the number of actually
 94        used halfedges, but corresponds to the size of the array allocated
 95        for halfedges."""
 96        return lib_py_gel.Manifold_no_allocated_halfedges(self.obj)
 97    def vertices(self):
 98        """ Returns an iterable containing all vertex indices"""
 99        verts = IntVector()
100        n = lib_py_gel.Manifold_vertices(self.obj, verts.obj)
101        return verts
102    def faces(self):
103        """ Returns an iterable containing all face indices"""
104        faces = IntVector()
105        n = lib_py_gel.Manifold_faces(self.obj, faces.obj)
106        return faces
107    def halfedges(self):
108        """ Returns an iterable containing all halfedge indices"""
109        hedges = IntVector()
110        n = lib_py_gel.Manifold_halfedges(self.obj, hedges.obj)
111        return hedges
112    def circulate_vertex(self, vid, mode='v'):
113        """ Circulate a vertex. Passed a vertex index, vid, and second argument,
114        mode='f', this function will return an iterable with all faces incident
115        on vid arranged in counter clockwise order. Similarly, if mode is 'h',
116        incident halfedges (outgoing) are returned, and for mode = 'v', all
117        neighboring vertices are returned. """
118        nbrs = IntVector()
119        n = lib_py_gel.Manifold_circulate_vertex(self.obj, vid, ct.c_char(mode.encode('ascii')), nbrs.obj)
120        return nbrs
121    def circulate_face(self, fid, mode='v'):
122        """ Circulate a face. Passed a face index, fid, and second argument,
123        mode='f', this function will return an iterable with all faces that
124        share an edge with fid (in counter clockwise order). If the argument is
125        mode='h', the halfedges themselves are returned. For mode='v', the
126        incident vertices of the face are returned. """
127        nbrs = IntVector()
128        n = lib_py_gel.Manifold_circulate_face(self.obj, fid, ct.c_char(mode.encode('ascii')), nbrs.obj)
129        return nbrs
130    def next_halfedge(self,hid):
131        """ Returns next halfedge to hid. """
132        return lib_py_gel.Walker_next_halfedge(self.obj, hid)
133    def prev_halfedge(self,hid):
134        """ Returns previous halfedge to hid. """
135        return lib_py_gel.Walker_prev_halfedge(self.obj, hid)
136    def opposite_halfedge(self,hid):
137        """ Returns opposite halfedge to hid. """
138        return lib_py_gel.Walker_opposite_halfedge(self.obj, hid)
139    def incident_face(self,hid):
140        """ Returns face corresponding to hid. """
141        return lib_py_gel.Walker_incident_face(self.obj, hid)
142    def incident_vertex(self,hid):
143        """ Returns vertex corresponding to (or pointed to by) hid. """
144        return lib_py_gel.Walker_incident_vertex(self.obj, hid)
145    def remove_vertex(self,vid):
146        """ Remove vertex vid from the Manifold. This function merges all faces
147        around the vertex into one and then removes this resulting face. """
148        return lib_py_gel.Manifold_remove_vertex(self.obj, vid)
149    def remove_face(self,fid):
150        """ Removes a face, fid, from the Manifold. If it is an interior face it is
151        simply replaced by an invalid index. If the face contains boundary
152        edges, these are removed. Situations may arise where the mesh is no
153        longer manifold because the situation at a boundary vertex is not
154        homeomorphic to a half disk. This, we can probably ignore since from the
155        data structure point of view it is not really a problem that a vertex is
156        incident on two holes - a hole can be seen as a special type of face.
157        The function returns false if the index of the face is not valid,
158        otherwise the function must complete. """
159        return lib_py_gel.Manifold_remove_face(self.obj, fid)
160    def remove_edge(self,hid):
161        """ Remove an edge, hid, from the Manifold. This function will remove the
162        faces on either side and the edge itself in the process. Thus, it is a
163        simple application of remove_face. """
164        return lib_py_gel.Manifold_remove_edge(self.obj, hid)
165    def vertex_in_use(self,vid):
166        """ check if vertex, vid, is in use. This function returns true if the id corresponds
167        to a vertex that is currently in the mesh and false otherwise. vid could
168        be outside the range of used ids and it could also correspond to a vertex
169        which is not active. The function returns false in both cases. """
170        return lib_py_gel.Manifold_vertex_in_use(self.obj, vid)
171    def face_in_use(self,fid):
172        """ check if face, fid, is in use. This function returns true if the id corresponds
173        to a face that is currently in the mesh and false otherwise. fid could
174        be outside the range of used ids and it could also correspond to a face
175        which is not active. The function returns false in both cases. """
176        return lib_py_gel.Manifold_face_in_use(self.obj, fid)
177    def halfedge_in_use(self,hid):
178        """ check if halfedge hid is in use. This function returns true if the id corresponds
179        to a halfedge that is currently in the mesh and false otherwise. hid could
180        be outside the range of used ids and it could also correspond to a halfedge
181        which is not active. The function returns false in both cases. """
182        return lib_py_gel.Manifold_halfedge_in_use(self.obj, hid)
183    def flip_edge(self,hid):
184        """ Flip the edge, hid, separating two faces. The function first verifies that
185        the edge is flippable. This entails making sure that all of the
186        following are true.
187        1. adjacent faces are triangles.
188        2. neither end point has valency three or less.
189        3. the vertices that will be connected are not already.
190        If the tests are passed, the flip is performed and the function
191        returns True. Otherwise False."""
192        return lib_py_gel.Manifold_flip_edge(self.obj,hid)
193    def collapse_edge(self,hid, avg_vertices=False):
194        """ Collapse an edge hid.
195        Before collapsing hid, a number of tests are made:
196        ---
197        1.  For the two vertices adjacent to the edge, we generate a list of all their neighbouring vertices.
198        We then generate a  list of the vertices that occur in both these lists.
199        That is, we find all vertices connected by edges to both endpoints of the edge and store these in a list.
200        2.  For both faces incident on the edge, check whether they are triangular.
201        If this is the case, the face will be removed, and it is ok that the the third vertex is connected to both endpoints.
202        Thus the third vertex in such a face is removed from the list generated in 1.
203        3.  If the list is now empty, all is well.
204        Otherwise, there would be a vertex in the new mesh with two edges connecting it to the same vertex. Return false.
205        4.  TETRAHEDRON TEST:
206        If the valency of both vertices is three, and the incident faces are triangles, we also disallow the operation.
207        Reason: A vertex valency of two and two triangles incident on the adjacent vertices makes the construction collapse.
208        5.  VALENCY 4 TEST:
209        If a triangle is adjacent to the edge being collapsed, it disappears.
210        This means the valency of the remaining edge vertex is decreased by one.
211        A valency two vertex reduced to a valency one vertex is considered illegal.
212        6.  PREVENT MERGING HOLES:
213        Collapsing an edge with boundary endpoints and valid faces results in the creation where two holes meet.
214        A non manifold situation. We could relax this...
215        7. New test: if the same face is in the one-ring of both vertices but not adjacent to the common edge,
216        then the result of a collapse would be a one ring where the same face occurs twice. This is disallowed as the resulting
217        face would be non-simple.
218        If the tests are passed, the collapse is performed and the function
219        returns True. Otherwise False."""
220        return lib_py_gel.Manifold_collapse_edge(self.obj, hid, avg_vertices)
221    def split_face_by_edge(self,fid,v0,v1):
222        """   Split a face. The face, fid, is split by creating an edge with
223        endpoints v0 and v1 (the next two arguments). The vertices of the old
224        face between v0 and v1 (in counter clockwise order) continue to belong
225        to fid. The vertices between v1 and v0 belong to the new face. A handle to
226        the new face is returned. """
227        return lib_py_gel.Manifold_split_face_by_edge(self.obj, fid, v0, v1)
228    def split_face_by_vertex(self,fid):
229        """   Split a polygon, fid, by inserting a vertex at the barycenter. This
230        function is less likely to create flipped triangles than the
231        split_face_triangulate function. On the other hand, it introduces more
232        vertices and probably makes the triangles more acute. The vertex id of the
233        inserted vertex is returned. """
234        return lib_py_gel.Manifold_split_face_by_vertex(self.obj,fid)
235    def split_edge(self,hid):
236        """   Insert a new vertex on halfedge hid. The new halfedge is insterted
237        as the previous edge to hid. The vertex id of the inserted vertex is returned. """
238        return lib_py_gel.Manifold_split_edge(self.obj,hid)
239    def stitch_boundary_edges(self,h0,h1):
240        """   Stitch two halfedges. Two boundary halfedges, h0 and h1, can be stitched
241        together. This can be used to build a complex mesh from a bunch of
242        simple faces. """
243        return lib_py_gel.Manifold_stitch_boundary_edges(self.obj, h0, h1)
244    def merge_faces(self,hid):
245        """   Merges two faces into a single polygon. The merged faces are those shared
246        by the edge for which hid is one of the two corresponding halfedges. This function returns
247        true if the merging was possible and false otherwise. Currently merge
248        only fails if the mesh is already illegal. Thus it should, in fact,
249        never fail. """
250        if self.is_halfedge_at_boundary(hid):
251            return False
252        fid = self.incident_face(hid)
253        return lib_py_gel.Manifold_merge_faces(self.obj, fid, hid)
254    def close_hole(self,hid):
255        """ Close hole given by hid (i.e. the face referenced by hid). Returns
256        index of the created face or the face that was already there if, in
257        fact, hid was not next to a hole. """
258        return lib_py_gel.Manifold_close_hole(self.obj, hid)
259    def cleanup(self):
260        """ Remove unused items from Mesh. This function remaps all vertices, halfedges
261        and faces such that the arrays do not contain any holes left by unused mesh
262        entities. It is a good idea to call this function when a mesh has been simplified
263        or changed in other ways such that mesh entities have been removed. However, note
264        that it invalidates any attributes that you might have stored in auxilliary arrays."""
265        lib_py_gel.Manifold_cleanup(self.obj)
266    def is_halfedge_at_boundary(self, hid):
267        """ Returns True if hid is a boundary halfedge, i.e. face on either
268        side is invalid. """
269        return lib_py_gel.is_halfedge_at_boundary(self.obj, hid)
270    def is_vertex_at_boundary(self, vid):
271        """ Returns True if vid lies on a boundary. """
272        return lib_py_gel.is_vertex_at_boundary(self.obj, vid)
273    def edge_length(self, hid):
274        """ Returns length of edge given by halfedge hid which is passed as argument. """
275        return lib_py_gel.length(self.obj, hid)
276    def valency(self,vid):
277        """ Returns valency of vid, i.e. number of incident edges."""
278        return lib_py_gel.valency(self.obj,vid)
279    def vertex_normal(self, vid):
280        """ Returns vertex normal (angle weighted) of vertex given by vid """
281        n = ndarray(3,dtype=np.float64)
282        lib_py_gel.vertex_normal(self.obj, vid, n)
283        return n
284    def connected(self, v0, v1):
285        """ Returns true if the two argument vertices, v0 and v1, are in each other's one-rings."""
286        return lib_py_gel.connected(self.obj,v0,v1)
287    def no_edges(self, fid):
288        """ Compute the number of edges of a face fid """
289        return lib_py_gel.no_edges(self.obj, fid)
290    def face_normal(self, fid):
291        """ Compute the normal of a face fid. If the face is not a triangle,
292        the normal is not defined, but computed using the first three
293        vertices of the face. """
294        n = ndarray(3, dtype=np.float64)
295        lib_py_gel.face_normal(self.obj, fid, n)
296        return n
297    def area(self, fid):
298        """ Returns the area of a face fid. """
299        return lib_py_gel.area(self.obj, fid)
300    def perimeter(self, fid):
301        """ Returns the perimeter of a face fid. """
302        return lib_py_gel.perimeter(self.obj, fid)
303    def centre(self, fid):
304        """ Returns the centre of a face. """
305        c = ndarray(3, dtype=np.float64)
306        lib_py_gel.centre(self.obj, fid, c)
307        return c
308
309def valid(m):
310    """This function performs a series of tests to check that this
311    is a valid manifold. This function is not rigorously constructed but seems
312    to catch all problems so far. The function returns true if the mesh is valid
313    and false otherwise. """
314    return lib_py_gel.valid(m.obj)
315
316def closed(m):
317    """ Returns true if m is closed, i.e. has no boundary."""
318    return lib_py_gel.closed(m.obj)
319
320def bbox(m):
321    """ Returns the min and max corners of the bounding box of Manifold m. """
322    pmin = ndarray(3,dtype=np.float64)
323    pmax = ndarray(3,dtype=np.float64)
324    lib_py_gel.bbox(m.obj, pmin, pmax)
325    return pmin, pmax
326
327def bsphere(m):
328    """ Calculate the bounding sphere of the manifold m.
329    Returns centre,radius """
330    c = ndarray(3,dtype=np.float64)
331    r = (ct.c_double)()
332    lib_py_gel.bsphere(m.obj, c, ct.byref(r))
333    return (c,r)
334
335def stitch(m, rad=1e-30):
336    """ Stitch together edges of m whose endpoints coincide geometrically. This
337    function allows you to create a mesh as a bunch of faces and then stitch
338    these together to form a coherent whole. What this function adds is a
339    spatial data structure to find out which vertices coincide. The return value
340    is the number of edges that could not be stitched. Often this is because it
341    would introduce a non-manifold situation."""
342    return lib_py_gel.stitch_mesh(m.obj,rad)
343
344def obj_save(fn, m):
345    """ Save Manifold m to Wavefront obj file. """
346    s = ct.c_char_p(fn.encode('utf-8'))
347    lib_py_gel.obj_save(s, m.obj)
348
349def off_save(fn, m):
350    """ Save Manifold m to OFF file. """
351    s = ct.c_char_p(fn.encode('utf-8'))
352    lib_py_gel.off_save(s, m.obj)
353
354def x3d_save(fn, m):
355    """ Save Manifold m to X3D file. """
356    s = ct.c_char_p(fn.encode('utf-8'))
357    lib_py_gel.x3d_save(s, m.obj)
358
359def obj_load(fn):
360    """ Load and return Manifold from Wavefront obj file.
361    Returns None if loading failed. """
362    m = Manifold()
363    s = ct.c_char_p(fn.encode('utf-8'))
364    if lib_py_gel.obj_load(s, m.obj):
365        return m
366    return None
367
368def off_load(fn):
369    """ Load and return Manifold from OFF file.
370    Returns None if loading failed."""
371    m = Manifold()
372    s = ct.c_char_p(fn.encode('utf-8'))
373    if lib_py_gel.off_load(s, m.obj):
374        return m
375    return None
376
377def ply_load(fn):
378    """ Load and return Manifold from Stanford PLY file.
379    Returns None if loading failed. """
380    m = Manifold()
381    s = ct.c_char_p(fn.encode('utf-8'))
382    if lib_py_gel.ply_load(s, m.obj):
383        return m
384    return None
385
386def x3d_load(fn):
387    """ Load and return Manifold from X3D file.
388    Returns None if loading failed."""
389    m = Manifold()
390    s = ct.c_char_p(fn.encode('utf-8'))
391    if lib_py_gel.x3d_load(s, m.obj):
392        return m
393    return None
394
395from os.path import splitext
396def load(fn):
397    """ Load a Manifold from an X3D/OBJ/OFF/PLY file. Return the
398    loaded Manifold. Returns None if loading failed."""
399    name, extension = splitext(fn)
400    if extension.lower() == ".x3d":
401        return x3d_load(fn)
402    if extension.lower() == ".obj":
403        return obj_load(fn)
404    if extension.lower() == ".off":
405        return off_load(fn)
406    if extension.lower() == ".ply":
407        return ply_load(fn)
408    return None
409
410def save(fn, m):
411    """ Save a Manifold, m, to an X3D/OBJ/OFF file. """
412    name, extension = splitext(fn)
413    if extension.lower() == ".x3d":
414        x3d_save(fn, m)
415    elif extension.lower() == ".obj":
416        obj_save(fn, m)
417    elif extension.lower() == ".off":
418        off_save(fn, m)
419
420
421def remove_caps(m, thresh=2.9):
422    """ Remove caps from a manifold, m, consisting of only triangles. A cap is a
423    triangle with two very small angles and an angle close to pi, however a cap
424    does not necessarily have a very short edge. Set the ang_thresh to a value
425    close to pi. The closer to pi the _less_ sensitive the cap removal. A cap is
426    removed by flipping the (long) edge E opposite to the vertex V with the
427    angle close to pi. However, the function is more complex. Read code and
428    document more carefully !!! """
429    lib_py_gel.remove_caps(m.obj,thresh)
430
431def remove_needles(m, thresh=0.05, average_positions=False):
432    """  Remove needles from a manifold, m, consisting of only triangles. A needle
433    is a triangle with a single very short edge. It is moved by collapsing the
434    short edge. The thresh parameter sets the length threshold (in terms of the average edge length
435    in the mesh). If average_positions is true then the collapsed vertex is placed at the average position of the end points."""
436    abs_thresh = thresh * average_edge_length(m)
437    lib_py_gel.remove_needles(m.obj,abs_thresh, average_positions)
438
439def close_holes(m, max_size=100):
440    """  This function replaces holes in m by faces. It is really a simple function
441    that just finds all loops of edges next to missing faces. """
442    lib_py_gel.close_holes(m.obj, max_size)
443
444def flip_orientation(m):
445    """  Flip the orientation of a mesh, m. After calling this function, normals
446    will point the other way and clockwise becomes counter clockwise """
447    lib_py_gel.flip_orientation(m.obj)
448
449def merge_coincident_boundary_vertices(m, rad = 1.0e-30):
450    """  Merge vertices of m that are boundary vertices and coincident.
451        However, if one belongs to the other's one ring or the one
452        rings share a vertex, they will not be merged. """
453    lib_py_gel.merge_coincident_boundary_vertices(m.obj, rad)
454
455def minimize_curvature(m,anneal=False):
456    """ Minimizes mean curvature of m by flipping edges. Hence, no vertices are moved.
457     This is really the same as dihedral angle minimization, except that we weight by edge length. """
458    lib_py_gel.minimize_curvature(m.obj, anneal)
459
460def minimize_dihedral_angle(m,max_iter=10000, anneal=False, alpha=False, gamma=4.0):
461    """ Minimizes dihedral angles in m by flipping edges.
462        Arguments:
463        max_iter is the maximum number of iterations for simulated annealing.
464        anneal tells us the code whether to apply simulated annealing
465        alpha=False means that we use the cosine of angles rather than true angles (faster)
466        gamma is the power to which the angles are raised."""
467    lib_py_gel.minimize_dihedral_angle(m.obj, max_iter, anneal,alpha,ct.c_double(gamma))
468
469
470def maximize_min_angle(m,dihedral_thresh=0.95,anneal=False):
471    """ Maximizes the minimum angle of triangles by flipping edges of m. Makes the mesh more Delaunay."""
472    lib_py_gel.maximize_min_angle(m.obj,dihedral_thresh,anneal)
473
474def optimize_valency(m,anneal=False):
475    """ Tries to achieve valence 6 internally and 4 along edges by flipping edges of m. """
476    lib_py_gel.optimize_valency(m.obj, anneal)
477
478def randomize_mesh(m,max_iter=1):
479    """  Make random flips in m. Useful for generating synthetic test cases. """
480    lib_py_gel.randomize_mesh(m.obj, max_iter)
481
482def quadric_simplify(m,keep_fraction,singular_thresh=1e-4,error_thresh=1):
483    """ Garland Heckbert simplification of mesh m. keep_fraction is the fraction of vertices
484    to retain. The singular_thresh determines how subtle features are preserved. For values
485    close to 1 the surface is treated as smooth even in the presence of sharp edges of low
486    dihedral angle (angle between normals). Close to zero, the method preserves even subtle
487    sharp features better. The error_thresh is the value of the QEM error at which
488    simplification stops. It is relative to the bounding box size. The default value is 1
489    meaning that simplification continues until the model has been simplified to a number of
490    vertices approximately equal to keep_fraction times the original number of vertices."""
491    lib_py_gel.quadric_simplify(m.obj, keep_fraction, singular_thresh,error_thresh)
492
493def average_edge_length(m):
494    """ Returns the average edge length of mesh m. """
495    return lib_py_gel.average_edge_length(m.obj)
496
497def median_edge_length(m):
498    """ Returns the median edge length of m"""
499    return lib_py_gel.median_edge_length(m.obj)
500
501def refine_edges(m,threshold):
502    """ Split all edges in m which are longer
503    than the threshold (second arg) length. A split edge
504    results in a new vertex of valence two."""
505    return lib_py_gel.refine_edges(m.obj, threshold)
506
507def cc_split(m):
508    """ Perform a Catmull-Clark split on m, i.e. a split where each face is divided
509    into new quadrilateral faces formed by connecting a corner with a point on
510    each incident edge and a point at the centre of the face."""
511    lib_py_gel.cc_split(m.obj)
512
513def loop_split(m):
514    """ Perform a loop split on m where each edge is divided into two segments, and
515    four new triangles are created for each original triangle. """
516    lib_py_gel.loop_split(m.obj)
517
518def root3_subdivide(m):
519    """ Leif Kobbelt's subdivision scheme applied to m. A vertex is placed in the
520    center of each face and all old edges are flipped. """
521    lib_py_gel.root3_subdivide(m.obj)
522
523def rootCC_subdivide(m):
524    """ This subdivision scheme creates a vertex inside each original (quad) face of m,
525    producing four triangles. Triangles sharing an old edge are then merged.
526    Two steps produce something similar to Catmull-Clark. """
527    lib_py_gel.rootCC_subdivide(m.obj)
528
529def butterfly_subdivide(m):
530    """ Butterfly subidiviosn on m. An interpolatory scheme. Creates the same connectivity as Loop. """
531    lib_py_gel.butterfly_subdivide(m.obj)
532
533def cc_smooth(m):
534    """ If called after cc_split, this function completes a step of Catmull-Clark
535    subdivision of m."""
536    lib_py_gel.cc_smooth(m.obj)
537
538def cc_subdivide(m):
539    """ Perform a full Catmull-Clark subdivision step on mesh m. """
540    lib_py_gel.cc_split(m.obj)
541    lib_py_gel.cc_smooth(m.obj)
542
543def loop_subdivide(m):
544    """ Perform a full Loop subdivision step on mesh m. """
545    lib_py_gel.loop_split(m.obj)
546    lib_py_gel.loop_smooth(m.obj)   
547
548def volume_preserving_cc_smooth(m, iter):
549    """ This function does the same type of smoothing as in Catmull-Clark
550    subdivision, but to preserve volume it actually performs two steps, and the
551    second step is negative as in Taubin smoothing."""
552    lib_py_gel.volume_preserving_cc_smooth(m.obj, iter)
553
554def regularize_quads(m, w=0.5, shrink=0.0, iter=1):
555    """ This function smooths a quad mesh by regularizing quads. Essentially,
556    regularization just makes them more rectangular. """
557    lib_py_gel.regularize_quads(m.obj, w, shrink, iter)
558
559
560def loop_smooth(m):
561    """ If called after Loop split, this function completes a step of Loop
562    subdivision of m. """
563    lib_py_gel.loop_smooth(m.obj)
564    
565def taubin_smooth(m, iter=1):
566    """ This function performs Taubin smoothing on the mesh m for iter number
567    of iterations. """
568    lib_py_gel.taubin_smooth(m.obj, iter)
569
570def laplacian_smooth(m, w=0.5, iter=1):
571    """ This function performs Laplacian smoothing on the mesh m for iter number
572    of iterations. w is the weight applied. """
573    lib_py_gel.laplacian_smooth(m.obj, w, iter)
574
575def anisotropic_smooth(m, sharpness=0.5, iter=1):
576    """ This function performs anisotropic smoothing on the mesh m for iter number
577    of iterations. A bilateral filtering controlled by sharpness is performed on 
578    the face normals followed by a rotation of the faces to match the new normals.
579    The updated vertex positions are the average positions of the corners of the 
580    rotated faces. For sharpness==0 the new normal is simply the area weighted 
581    average of the normals of incident faces. For sharpness>0 the weight of the
582    neighbouring face normals is a Gaussian function of the angle between the
583    face normals. The greater the sharpness, the more the smoothing is anisotropic."""
584    lib_py_gel.anisotropic_smooth(m.obj, sharpness, iter)
585    
586def volumetric_isocontour(data, bbox_min = None, bbox_max = None,
587                          tau=0.0,
588                          make_triangles=True,
589                          high_is_inside=True,
590                          dual_connectivity=False):
591    """ Creates a polygonal mesh from volumetric data by isocontouring. The dimensions
592    are given by dims, bbox_min (defaults to [0,0,0] ) and bbox_max (defaults to dims) are
593    the corners of the bounding box in R^3 that corresponds to the volumetric grid, tau is
594    the iso value (defaults to 0). If make_triangles is True (default), we turn the quads
595    into triangles. Finally, high_is_inside=True (default) means that values greater than
596    tau are interior and smaller values are exterior. If dual_connectivity is False (default)
597    the function produces marching cubes connectivity and if True it produces dual contouring
598    connectivity. MC connectivity tends to produce less nice triangle shapes but since the
599    vertices always lie on edges, the geometry is arguably better defined for MC. """
600    m = Manifold()
601    dims = data.shape
602    if bbox_min is None:
603        bbox_min = (0,0,0)
604    if bbox_max is None:
605        bbox_max = dims
606    data_float = np.asarray(data, dtype=ct.c_float, order='F')
607    bbox_min_d = np.asarray(bbox_min, dtype=np.float64, order='C')
608    bbox_max_d = np.asarray(bbox_max, dtype=np.float64, order='C')
609    lib_py_gel.volumetric_isocontour(m.obj, dims[0], dims[1], dims[2],
610                                     data_float, bbox_min_d, bbox_max_d, tau,
611                                     make_triangles, high_is_inside, dual_connectivity)
612    return m
613
614def triangulate(m, clip_ear=True):
615    """ Turn a general polygonal mesh, m, into a triangle mesh by repeatedly
616        splitting a polygon into smaller polygons. """
617    if clip_ear:
618        lib_py_gel.ear_clip_triangulate(m.obj)
619    else:
620        lib_py_gel.shortest_edge_triangulate(m.obj)
621
622def skeleton_to_feq(g, node_radii = None, symmetrize=True):
623    """ Turn a skeleton graph g into a Face Extrusion Quad Mesh m with given node_radii for each graph node.
624    If symmetrize is True (default) the graph is made symmetrical. If node_radii are supplied then they
625    are used in the reconstruction. Otherwise, the radii are obtained from the skeleton. They are stored in 
626    the green channel of the vertex color during skeletonization, so for a skeletonized shape that is how the
627    radius of each node is obtained. This is a questionable design decision and will probably change 
628    in the future. """
629    m = Manifold()
630    if node_radii is None:
631        node_radii = [0.0] * len(g.nodes())
632        use_graph_radii = True
633    else:
634        use_graph_radii = False
635        if isinstance(node_radii, (int, float)):
636            if node_radii <= 0.0:
637                node_radii = 0.25 * g.average_edge_length()
638            node_radii = [node_radii] * len(g.nodes())
639
640    node_rs_flat = np.asarray(node_radii, dtype=np.float64)
641    lib_py_gel.graph_to_feq(g.obj , m.obj, node_rs_flat, symmetrize, use_graph_radii)
642    return m
643
644# def non_rigid_registration(m, ref_mesh):
645#     """ Perform non-rigid registration of m to ref_mesh. """
646#     lib_py_gel.non_rigid_registration(m.obj, ref_mesh.obj)
647
648def laplacian_matrix(m):
649    num_verts = m.no_allocated_vertices()
650    laplacian = np.full((num_verts,num_verts), 0.0)
651    for i in m.vertices():
652        neighbors = m.circulate_vertex(i)
653        deg = len(neighbors)
654        laplacian[i][i] = 1.0
655        for v in neighbors:
656            laplacian[i][v] = -1/deg
657    return csc_matrix(laplacian)
658
659
660def inv_correspondence_leqs(m, ref_mesh):
661    m_pos = m.positions()
662    ref_pos = ref_mesh.positions()
663
664    m_tree = spatial.I3DTree()
665    for v in m.vertices():
666        m_tree.insert(m_pos[v],v)
667    m_tree.build()
668
669    m_target_pos = np.zeros(m.positions().shape)
670    m_cnt = np.zeros(m.no_allocated_vertices())
671    for r_id in ref_mesh.vertices():
672        query_pt = ref_pos[r_id]
673        closest_pt_obj = m_tree.closest_point(query_pt,1e32)
674        if closest_pt_obj is not None:
675            key,m_id = closest_pt_obj
676            r_norm = ref_mesh.vertex_normal(r_id)
677            m_norm = m.vertex_normal(m_id)
678            dot_prod = m_norm @ r_norm
679            if dot_prod > 0.5:
680                v = query_pt - m_pos[m_id]
681                wgt = dot_prod - 0.5
682                m_target_pos[m_id] += wgt*query_pt
683                m_cnt[m_id] += wgt
684                
685    N = m.no_allocated_vertices()
686    A_list = []
687    b_list = []
688    for vid in m.vertices():
689        if m_cnt[vid] > 0.0:
690            row_a = np.zeros(N)
691            row_a[vid] = 1.0
692            A_list.append(row_a)
693            b_list.append(m_target_pos[vid]/m_cnt[vid])
694    return csc_matrix(np.array(A_list)), np.array(b_list)
695
696def fit_mesh_to_ref(m, ref_mesh, iter = 10, dist_wt = 1.0, lap_wt = 0.3):
697    """ Fits a skeletal mesh m to a reference mesh ref_mesh. """
698
699    v_pos = m.positions()
700    for i in range(iter):
701        Ai, bi = inv_correspondence_leqs(m, ref_mesh)
702        lap_matrix = laplacian_matrix(m)
703        lap_b = lap_matrix @ v_pos
704        final_A = vstack([lap_wt*lap_matrix, dist_wt*Ai])
705        final_b = np.vstack([lap_wt*lap_b, dist_wt*bi])
706        opt_x, _, _, _ = lsqr(final_A, final_b[:,0])[:4]
707        opt_y, _, _, _ = lsqr(final_A, final_b[:,1])[:4]
708        opt_z, _, _, _ = lsqr(final_A, final_b[:,2])[:4]
709        v_pos[:,:] = np.stack([opt_x, opt_y, opt_z], axis=1)
710        regularize_quads(m, w=0.5, shrink=0.8, iter=3)
711    return m
712
713class MeshDistance:
714    """ This class allows you to compute the distance from any point in space to
715    a Manifold (which must be triangulated). The constructor creates an instance
716    based on a specific mesh, and the signed_distance function computes the actual distance. """
717    def __init__(self,m):
718        self.obj = lib_py_gel.MeshDistance_new(m.obj)
719    def __del__(self):
720        lib_py_gel.MeshDistance_delete(self.obj)
721    def signed_distance(self,pts,upper=1e30):
722        """ Compute the signed distance from each point in pts to the mesh stored in
723        this class instance. pts should be convertible to a length N>=1 array of 3D
724        points. The function returns an array of N distance values with a single distance
725        for each point. The distance corresponding to a point is positive if the point
726        is outside and negative if inside. The upper parameter can be used to threshold
727        how far away the distance is of interest. """
728        p = np.asarray(pts, dtype=ct.c_float, order='C')
729        ndim = len(p.shape)
730        if ndim==1:
731            n = p.shape[0]//3
732        elif ndim==2:
733            n = p.shape[0]
734        else:
735            raise Exception("you must pass signed_distance pts as a 1D array or a 2D array of dim nx3")
736        
737        d = np.ndarray(n, dtype=ct.c_float)
738        lib_py_gel.MeshDistance_signed_distance(self.obj, n, p, d, upper)
739        return d
740    def ray_inside_test(self,pts,no_rays=3):
741        """Check whether each point in pts is inside or outside the stored mesh by
742        casting rays. pts should be convertible to a length N>=1 array of 3D points.
743        Effectively, this is the sign of the distance. In some cases casting (multiple)
744        ray is more robust than using the sign computed locally. Returns an array of
745        N integers which are either 1 or 0 depending on whether the corresponding point
746        is inside (1) or outside (0). """
747        p = np.asarray(pts, dtype=ct.c_float, order='C')
748        ndim = len(p.shape)
749        if ndim==1:
750            n = p.shape[0]//3
751        elif ndim==2:
752            n = p.shape[0]
753        else:
754            raise Exception("you must pass signed_distance pts as a 1D array or a 2D array of dim nx3")
755        s = np.ndarray(n, dtype=ct.c_int)
756        lib_py_gel.MeshDistance_ray_inside_test(self.obj,n,p,s,no_rays)
757        return s
758    def intersect(self, p0, dir, _t=0):
759        """ Intersect the ray starting in p0 with direction, dir, with the stored mesh. Returns
760        the point of intersection if there is one, otherwise None. """
761        p0 = np.asarray(p0,dtype=ct.c_float)
762        dir = np.asarray(dir,dtype=ct.c_float)
763        t = ct.c_float(_t)
764        r = lib_py_gel.MeshDistance_ray_intersect(self.obj, p0, dir, ct.byref(t))
765        if r:
766            return t.value, p0, dir
767        return None