skeletor.pre

 37def fix_mesh(mesh, remote_infinite=True, merge_duplicate_verts=True,
 38             remove_degenerate_faces=True, remove_unreferenced_verts=True,
 39             drop_winglets=True, fix_normals=False, remove_disconnected=False,
 40             inplace=False):
 41    """Try to fix some common problems with mesh.
 42
 43     1. Remove infinite values
 44     2. Merge duplicate vertices
 45     3. Remove duplicate and degenerate faces
 46     4. Remove unreference vertices
 47     5. Drop winglets (faces that have only one adjacent face)
 48     5. Fix normals (Optional)
 49     6. Remove disconnected fragments (Optional)
 50
 51    Parameters
 52    ----------
 53    mesh :                  mesh-like object
 54                            Mesh to fix. Must have `.vertices` and `.faces`
 55                            properties.
 56    remove_disconnected :   False | int
 57                            If a number is given, will iterate over the mesh's
 58                            connected components and remove those consisting of
 59                            less than the given number of vertices. For example,
 60                            ``remove_disconnected=5`` will drop parts of the
 61                            mesh that consist of five or less connected
 62                            vertices.
 63    inplace :               bool
 64                            If True, will perform fixes on the input mesh.
 65                            If False, will make a copy first. This is silently
 66                            ignored if `mesh` is not already a trimesh.
 67
 68    Returns
 69    -------
 70    fixed mesh :        trimesh.Trimesh
 71
 72    """
 73    if not isinstance(mesh, tm.Trimesh):
 74        mesh = make_trimesh(mesh, validate=False)
 75    elif not inplace:
 76        mesh = mesh.copy()
 77
 78    if remove_disconnected:
 79        to_drop = []
 80        for c in nx.connected_components(mesh.vertex_adjacency_graph):
 81            if len(c) <= remove_disconnected:
 82                to_drop += list(c)
 83
 84        # Remove dropped vertices
 85        remove = np.isin(np.arange(mesh.vertices.shape[0]), to_drop)
 86        mesh.update_vertices(~remove)
 87
 88    if remote_infinite:
 89        mesh.remove_infinite_values()
 90
 91    if merge_duplicate_verts:
 92        mesh.merge_vertices()
 93
 94    if remove_degenerate_faces:
 95        mesh.remove_duplicate_faces()
 96        mesh.remove_degenerate_faces()
 97
 98    if remove_unreferenced_verts:
 99        mesh.remove_unreferenced_vertices()
100
101    if drop_winglets:
102        mesh = remove_winglets(mesh)
103
104    # This should be done after all clean-up operations
105    if fix_normals:
106        mesh.fix_normals()
107
108    return mesh

192def simplify(mesh, ratio):
193    """Simplify mesh using Blender 3D.
194
195    Uses Blender's "decimate" modifier in "collapse" mode.
196
197    Parameters
198    ----------
199    mesh :  trimesh.Trimesh
200            Mesh to simplify.
201    ratio : float
202            Factor to which to reduce faces. For example, a ratio of 0.5 will
203            reduce the number of faces to 50%.
204
205    Returns
206    -------
207    trimesh.Trimesh
208            Simplified mesh.
209
210    """
211    if not tm.interfaces.blender.exists:
212        raise ImportError('No Blender available (executable not found).')
213    _blender_executable = tm.interfaces.blender._blender_executable
214
215    assert ratio < 1 and ratio > 0, 'ratio must be between 0 and 1'
216
217    # We need to import here to avoid circular imports
218    from ..utilities import make_trimesh
219    mesh = make_trimesh(mesh, validate=False)
220    assert isinstance(mesh, tm.Trimesh)
221
222    # Load the template
223    temp_name = 'blender_decimate.py.template'
224    if temp_name in _cache:
225        template = _cache[temp_name]
226    else:
227        with open(os.path.join(_pwd, 'templates', temp_name), 'r') as f:
228            template = f.read()
229        _cache[temp_name] = template
230
231    # Replace placeholder with actual ratio
232    script = template.replace('$RATIO', str(ratio))
233
234    # Let trimesh's MeshScript take care of exectution and clean-up
235    with tm.interfaces.generic.MeshScript(meshes=[mesh],
236                                          script=script,
237                                          debug=False) as blend:
238        result = blend.run(_blender_executable
239                           + ' --background --python $SCRIPT')
240
241    # Blender apparently returns actively incorrect face normals
242    result.face_normals = None
243
244    return result

247def remesh(mesh, voxel_size=50, adaptivity=5):
248    """Remesh mesh using Blender 3D.
249
250    Uses Blender's "remesh" modifier in "voxel" mode.
251
252    Parameters
253    ----------
254    mesh :  trimesh.Trimesh
255            Mesh to remesh.
256    voxel_size : float
257            Size of individual voxels (edge length).
258    adaptivity: float
259            Reduces final face count where detail is not important.
260
261    Returns
262    -------
263    trimesh.Trimesh
264            Remeshed mesh.
265
266    """
267    if not tm.interfaces.blender.exists:
268        raise ImportError('No Blender available (executable not found).')
269    _blender_executable = tm.interfaces.blender._blender_executable
270
271    assert voxel_size > 0, '`voxel_size` must be a positive number'
272    assert adaptivity > 0, '`adaptivity` must be a positive number'
273
274    # We need to import here to avoid circular imports
275    from ..utilities import make_trimesh
276    mesh = make_trimesh(mesh, validate=False)
277    assert isinstance(mesh, tm.Trimesh)
278
279    # Load the template
280    temp_name = 'blender_remesh.py.template'
281    if temp_name in _cache:
282        template = _cache[temp_name]
283    else:
284        with open(os.path.join(_pwd, 'templates', temp_name), 'r') as f:
285            template = f.read()
286        _cache[temp_name] = template
287
288    # Replace placeholder with actual ratio
289    script = template.replace('$VOXEL_SIZE', str(voxel_size)
290                              ).replace('$ADAPTIVITY', str(adaptivity))
291
292    # Let trimesh's MeshScript take care of exectution and clean-up
293    with tm.interfaces.generic.MeshScript(meshes=[mesh],
294                                          script=script,
295                                          debug=False) as blend:
296        result = blend.run(_blender_executable
297                           + ' --background --python $SCRIPT')
298
299    # Blender apparently returns actively incorrect face normals
300    result.face_normals = None
301
302    return result

 37def contract(mesh, epsilon=1e-06, iter_lim=100, time_lim=None, precision=1e-07,
 38             SL=2, WH0=1, WL0='auto', operator='cotangent', progress=True,
 39             validate=True):
 40    """Contract mesh.
 41
 42    In a nutshell: this function contracts the mesh by applying rounds of
 43    _constraint_ Laplacian smoothing. This function can be fairly expensive
 44    and I highly recommend you play around with ``SL`` to contract the
 45    mesh in as few steps as possible. The contraction doesn't have to be perfect
 46    for good skeletonization results (<10%, i.e. `epsilon<0.1`).
 47
 48    Also: parameterization matters a lot! Default parameters will get you there
 49    but playing around with `SL` and `WH0` might speed things up by an order of
 50    magnitude.
 51
 52    Parameters
 53    ----------
 54    mesh :          mesh obj
 55                    The mesh to be contracted. Can be any object (e.g.
 56                    a trimesh.Trimesh) that has ``.vertices`` and ``.faces``
 57                    properties or a tuple ``(vertices, faces)`` or a dictionary
 58                    ``{'vertices': vertices, 'faces': faces}``.
 59                    Vertices and faces must be (N, 3) numpy arrays.
 60    epsilon :       float (0-1), optional
 61                    Target contraction rate as measured by the sum of all face
 62                    areas in the contracted versus the original mesh. Algorithm
 63                    will stop once mesh is contracted below this threshold.
 64                    Depending on your mesh (number of faces, shape) reaching a
 65                    strong contraction can be extremely costly with comparatively
 66                    little benefit for the subsequent skeletonization. Note that
 67                    the algorithm might stop short of this target if ``iter_lim``
 68                    or ``time_lim`` is reached first or if the sum of face areas
 69                    is increasing from one iteration to the next instead of
 70                    decreasing.
 71    iter_lim :      int (>1), optional
 72                    Maximum rounds of contractions.
 73    time_lim :      int, optional
 74                    Maximum run time in seconds. Note that this limit is not
 75                    checked during but after each round of contraction. Hence,
 76                    the actual total time will likely overshoot ``time_lim``.
 77    precision :     float, optional
 78                    Sets the precision for finding the least-square solution.
 79                    This is the main determinant for speed vs quality: lower
 80                    values will take (much) longer but will get you closer to an
 81                    optimally contracted mesh. Higher values will be faster but
 82                    the iterative contractions might stop early.
 83    SL :            float, optional
 84                    Factor by which the contraction matrix is multiplied for
 85                    each iteration. Higher values = quicker contraction, lower
 86                    values = more likely to get you an optimal contraction.
 87    WH0 :           float, optional
 88                    Initial weight factor for the attraction constraints.
 89                    The ratio of the initial weights ``WL0`` and ``WH0``
 90                    controls the smoothness and the degree of contraction of the
 91                    first iteration result, thus it determines the amount of
 92                    details retained in subsequent and final contracted meshes:
 93                    higher ``WH0`` = more details retained.
 94    WL0 :           "auto" | float
 95                    Initial weight factor for the contraction constraints. By
 96                    default ("auto"), this will be set to ``1e-3 * sqrt(A)``
 97                    with ``A`` being the average face area. This ensures that
 98                    contraction forces scale with the coarseness of the mesh.
 99    operator :      "cotangent" | "umbrella"
100                    Which Laplacian operator to use:
101
102                      - The "cotangent" operator (default) takes both topology
103                        and geometry of the mesh into account and is hence a
104                        better descriptor of the curvature flow. This is the
105                        operator used in the original paper.
106                      - The "umbrella" operator (aka "uniform weighting") uses
107                        only topological features of the mesh. This also makes
108                        it more robust against flaws in the mesh! Use it when
109                        the cotangent operator produces oddly contracted meshes.
110
111    progress :      bool
112                    Whether or not to show a progress bar.
113    validate :      bool
114                    If True, will try to fix potential issues with the mesh
115                    (e.g. infinite values, duplicate vertices, degenerate faces)
116                    before collapsing. Degenerate meshes can lead to effectively
117                    infinite runtime for this function!
118
119    Returns
120    -------
121    trimesh.Trimesh
122                    Contracted copy of original mesh. The final contraction rate
123                    is attached to the mesh as ``.epsilon`` property.
124
125    References
126    ----------
127    [1] Au OK, Tai CL, Chu HK, Cohen-Or D, Lee TY. Skeleton extraction by mesh
128        contraction. ACM Transactions on Graphics (TOG). 2008 Aug 1;27(3):44.
129
130    """
131    assert operator in ('cotangent', 'umbrella')
132    start = time.time()
133
134    # Force into trimesh
135    m = make_trimesh(mesh, validate=validate)
136    n = len(m.vertices)
137
138    # Initialize attraction weights
139    zeros = np.zeros((n, 3))
140    WH0_diag = np.zeros(n)
141    WH0_diag.fill(WH0)
142    WH0 = sp.sparse.spdiags(WH0_diag, 0, WH0_diag.size, WH0_diag.size)
143    # Make a copy but keep original values
144    WH = sp.sparse.dia_matrix(WH0)
145
146    # Initialize contraction weights
147    if WL0 == 'auto':
148        WL0 = 1e-03 * np.sqrt(averageFaceArea(m))
149        #WL0 = 1.0 / 10.0 * np.sqrt(averageFaceArea(m))
150    WL_diag = np.zeros(n)
151    WL_diag.fill(WL0)
152    WL = sp.sparse.spdiags(WL_diag, 0, WL_diag.size, WL_diag.size)
153
154    # Copy mesh
155    dm = m.copy()
156
157    area_ratios = [1.0]
158    originalRingAreas = getOneRingAreas(dm)
159    goodvertices = dm.vertices
160    bar_format = ("{l_bar}{bar}| [{elapsed}<{remaining}, "
161                  "{postfix[0]}/{postfix[1]}it, "
162                  "{rate_fmt}, epsilon {postfix[2]:.2g}")
163    with tqdm(total=100,
164              bar_format=bar_format,
165              disable=progress is False,
166              postfix=[1, iter_lim, 1]) as pbar:
167        for i in range(iter_lim):
168            # Get Laplace weights
169
170            if operator == 'cotangent':
171                L = laplacian_cotangent(dm, normalized=True)
172            else:
173                L = laplacian_umbrella(dm)
174            """
175            import robust_laplacian
176            L, M_ = robust_laplacian.mesh_laplacian(np.array(dm.vertices),
177                                                    np.array(dm.faces),
178                                                    mollify_factor=1e-3)
179            """
180            V = getMeshVPos(dm)
181            A = sp.sparse.vstack([WL.dot(L), WH])
182            b = np.vstack((zeros, WH.dot(V)))
183
184            cpts = np.zeros((n, 3))
185            for j in range(3):
186                """
187                # Solve A*x = b
188                # Note that we force scipy's lsqr() to use current vertex
189                # positions as start points - this speeds things up and
190                # without it we get suboptimal solutions that lead to early
191                # termination
192                cpts[:, j] = lsqr(A, b[:, j],
193                                  atol=precision, btol=precision,
194                                  damp=1,
195                                  x0=dm.vertices[:, j])[0]
196                """
197                # The solution below is recommended in scipy's lsqr docstring
198                # for when we have an initial estimate
199                # Gives use the same results as above but is slightly faster
200
201                # Initial estimate (i.e. our current positions)
202                x0 = dm.vertices[:, j]
203                # Compute residual vector
204                r0 = b[:, j] - A * x0
205                # Use LSQR to solve the system
206                dx = lsqr(A, r0,
207                          atol=precision, btol=precision,
208                          damp=1)[0]
209                # Add the correction dx to obtain a final solution
210                cpts[:, j] = x0 + dx
211
212            # Update mesh with new vertex position
213            dm.vertices = cpts
214
215            # Update iteration in progress bar
216            if progress:
217                pbar.postfix[0] = i + 1
218
219            # Break if face area has increased compared to the last iteration
220            new_eps = dm.area / m.area
221            if (new_eps > area_ratios[-1]):
222                dm.vertices = goodvertices
223                if progress:
224                    tqdm.write("Total face area increased from last iteration."
225                               f" Contraction stopped prematurely after {i} "
226                               f"iterations at epsilon {area_ratios[-1]:.2g}.")
227                break
228            area_ratios.append(new_eps)
229
230            # Update progress bar
231            if progress:
232                pbar.postfix[2] = area_ratios[-1]
233                prog = round((area_ratios[-2] - area_ratios[-1]) / (1 - epsilon) * 100)
234                pbar.update(min(prog, 100-pbar.n))
235
236            goodvertices = cpts
237
238            # Update contraction weights -> at each iteration the contraction
239            # forces increase to counteract the increased attraction forces
240            WL = sp.sparse.dia_matrix(WL.multiply(SL))
241
242            # Update attraction weights -> the smaller the one ring areas
243            # the higher the attraction forces
244            changeinarea = np.sqrt(originalRingAreas / getOneRingAreas(dm))
245            WH = sp.sparse.dia_matrix(WH0.multiply(changeinarea))
246
247            # Stop if we reached our target contraction rate
248            if (area_ratios[-1] <= epsilon):
249                break
250
251            # Stop if time limit is reached
252            if not isinstance(time_lim, (bool, type(None))):
253                if (time.time() - start) >= time_lim:
254                    break
255
256        # Keep track of final epsilon
257        dm.epsilon = area_ratios[-1]
258
259        return dm