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