Extract center curves from 3D meshes

I’m looking for a way to extract the skeleton of a mesh - the center curves of the solid mesh.
I have previously used ngon mesh skeleton component, but in the latest version this component is no longer there. (@Petras_Vestartas).

Are there any alternatives? Or suggestions for methods.
This article describes pretty well what I intend to do

Image from article:

Recently I added this method to my compas_wood plugin:

#! python3
# venv: timber_connections
# r: wood-nano==0.2.0, compas-wood==2.4.0, compas==2.8.1, wood-rui==0.2.2

import Rhino
import Rhino.Input
import Rhino.DocObjects
from Rhino.Geometry import Plane
from typing import *
from wood_rui import wood_rui_globals, process_input, add_skeleton
from compas_wood.binding_rhino import beam_skeleton, mesh_skeleton
import System


def callback(selection: dict[str, any], name: str):
    # Just axis without radii
    if not selection["simplify"]:
        for mesh in selection["meshes"]:
            polylines = mesh_skeleton(mesh)
            for polyline in polylines:
                Rhino.RhinoDoc.ActiveDoc.Objects.AddPolyline(polyline)
        return

    # Try to simplify an axis

    polylines = []
    distances = []
    meshes = []
    transforms = []

    for mesh in selection["meshes"]:
        polyline = None
        polylines = []

        polyline, distance = beam_skeleton(
            mesh,
            selection["divisions"],
            selection["nearest_neighbors"],
            selection["extend_end"],
        )
        polylines.append(polyline)
        distances.append(distance)
        meshes.append(mesh)

        # Add a small segement at the end of the polyline to provide orientation
        z_axis = polyline[polyline.Count - 1] - polyline[0]
        plane = Plane(polyline[polyline.Count - 1], z_axis)
        x_axis = plane.XAxis
        polyline.Add(polyline[polyline.Count - 1] + x_axis * 0.1)
        plane_xy = Plane.WorldXY
        plane = Plane(
            polyline[polyline.Count - 1],
            x_axis,
            Rhino.Geometry.Vector3d.CrossProduct(x_axis, z_axis),
        )
        xform = Rhino.Geometry.Transform.PlaneToPlane(plane, plane_xy)
        transforms.append(xform)

    add_skeleton(polylines, "mesh_skeleton", distances, meshes, transforms)


if __name__ == "__main__":
    input_dict = {
        "meshes": ([], list[Rhino.Geometry.Mesh]),
        "simplify": (False, bool),
        "divisions": (10, int),
        "nearest_neighbors": (10, int),
        "extend_end": (True, bool),
    }

    # Call the generalized input method with the dataset name and input dictionary
    process_input(input_dict, callback, hide_input=False)

This is method what you up to, since you dont care about straight timber beams made from scans. This method should work on your favorite editors rhino script editor, grasshopper or vscode:
from compas_wood.binding_rhino import mesh_skeleton

That is irrelated of the shape:

CGAL docs:

All the credits must be given to the smart people:

The initial implementation of this package is the result of the work of Xiang Gao during the 2013 season of the Google Summer of Code mentored by Andrea Tagliasacchi and Sébastien Loriot. It was finalized by Andreas Fabri and Sébastien Loriot.

At ibois long time ago I learnt how to wrap C++ methods to C# and Python since all the goodies are in the low level languages…

You are fantastic - absolutely fantastic. Thank you!

So I tried implementing this in a python component but the result is not easy to get as great as it was in ngon.

This is my code

#!/usr/bin/env python3
# venv: timber_connections
# r: wood-nano==0.2.0, compas-wood==2.4.0, compas==2.8.1, wood-rui==0.2.2


import System
import Rhino
import Rhino.Input
import Rhino.DocObjects
from Rhino.Geometry import Plane
from typing import *
from wood_rui import wood_rui_globals, process_input, add_skeleton
from compas_wood.binding_rhino import beam_skeleton, mesh_skeleton
import System

def compute_mesh_skeleton(input_meshes: list[Rhino.Geometry.Mesh],
                          simplify: bool = True,
                          divisions: int = 100,
                          nearest_neighbors: int = 10,
                          extend_end: bool = False
                         ) -> dict[str, any]:
    if not simplify:
        # Use mesh_skeleton when not simplifying: accumulate all returned polylines.
        polylines = []
        for mesh in input_meshes:
            for poly in mesh_skeleton(mesh):
                polylines.append(poly)
        return {"polylines": polylines}
    
    # When simplifying, use beam_skeleton and adjust each polyline.
    polylines = []
    distances = []
    meshes_out = []
    transforms = []
    
    for mesh in input_meshes:
        polyline, distance = beam_skeleton(mesh, divisions, nearest_neighbors, extend_end)
        polylines.append(polyline)
        distances.append(distance)
        meshes_out.append(mesh)
        
        # Add a small segment at the end to provide orientation.
        z_axis = polyline[polyline.Count - 1] - polyline[0]
        plane = Plane(polyline[polyline.Count - 1], z_axis)
        x_axis = plane.XAxis
        polyline.Add(polyline[polyline.Count - 1] + x_axis * 0.1)
        plane_xy = Plane.WorldXY
        plane = Plane(polyline[polyline.Count - 1],
                      x_axis,
                      Rhino.Geometry.Vector3d.CrossProduct(x_axis, z_axis))
        xform = Rhino.Geometry.Transform.PlaneToPlane(plane, plane_xy)
        transforms.append(xform)
    
    return {
        "polylines": polylines,
        "distances": distances,
        "meshes": meshes_out,
        "transforms": transforms
    }

skeleton = compute_mesh_skeleton([mesh])['polylines']

And here is the grasshopper script
MeshSkeleton.gh (957.5 KB)

Dense meshes:

Skeleton:

Dont use beam_skeleton but mesh_skeleton.

Try it and let me know.

polylines = mesh_skeleton(mesh)
            for polyline in polylines:
                Rhino.RhinoDoc.ActiveDoc.Objects.AddPolyline(polyline

Distances you can find by closest point to mesh search in rhinocommon.

I see. If I use simplify = false I get the mesh_skeleton function, but what it returns to me is this:

(I have used shrinkwrap to create a mesh)
MeshSkeleton.gh (962.4 KB)

Triangulate the mesh.

This is literally one line of code:

By the way, if you are using only 3d boxes, you could just compute principal-component-analysis from opennest. The normal of a plane is your axis:


THANK YOU!

Hi Petras,

First, thank you for sharing your extensive investment in time and expertise building these workflow automation tools. I’m looking to learn and be more of a computational designer, and I’m finding your forum and code contributions really helpful.

I just downloaded the MeshSkeleton example provided above and the python script 3 component is failing with the following error:

Error building script: Could not install packages due to an OSError: HTTPSConnectionPool(host=‘files.pythonhosted.org’, port=443): Max retries exceeded with url: /packages/8a/e0/c83016b086dbc68030d9617d9214411da232da26938eb9813f5ec2b4d394/wood_nano-0.2.0-cp39-cp39-win_amd64.whl.metadata (Caused by SSLError(SSLCertVerificationError(1, ‘[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed: self signed certificate in certificate chain (_ssl.c:1129)’)))

I suspect that given the way IT manages my work computer this may be a VPN related thing, but in reviewing the script itself I don’t see reference to wood-nano other than in the opening comments, and I haven’t studied the package dependencies at all yet (and I’m pretty new to scripting, so still finding my footing here…) so I’m not sure why this is reaching out to the internet for something.

Can you advise a workaround or help me understand what I need to do to get this running? Thanks!

I also implemented it here, can you try this:

just add in python r: compas, compas_cgal , cgal is pip installable so should work on rhino

hi @Petras_Vestartas

I could not find the Grasshopper Python component related

to row find Mesh Axis in your examples for the compas_wood plugin . Could you please share this component as a .gh file here?

@Petras_Vestartas :backhand_index_pointing_up:t3:

Have you tried to use this, mesh has to be closed and triangulated. Open grasshopper python plugin copy paste this code, and define input as a mesh:

#! python3

#venv: timber_connections

#r: wood-nano==0.2.0, compas-wood==2.4.0, compas==2.4.3, wood-rui==0.2.0

from compas_wood.binding_rhino import mesh_skeleton

polylines = mesh_skeleton(mesh)

Thanks @Petras_Vestartas

-give Error=``Error building script: Cannot install compas-wood==2.4.0 and compas==2.4.3 because these package versions have conflicting dependencies.```

and Error building script: Cannot install compas-wood==2.4.0 and compas==2.4.3 because these package versions have conflicting dependencies.

I have corrected that kind of code

#! python3
#venv: timber_connections
#r: compas==2.8.1, compas-wood==2.4.0, wood-nano==0.2.0, wood-rui==0.2.0

import subprocess
import sys

def check_version(pkg, target_version):
    try:
        import importlib.metadata as metadata
        version = metadata.version(pkg)
        return version == target_version
    except ImportError:
        return False
    except metadata.PackageNotFoundError:
        return False

def install_package(pkg):
    subprocess.check_call([sys.executable, "-m", "pip", "install", "--upgrade", pkg])

def uninstall_package(pkg):
    subprocess.call([sys.executable, "-m", "pip", "uninstall", "-y", pkg])

packages = {
    "compas": "2.8.1",
    "compas-wood": "2.4.0",
    "wood-nano": "0.2.0",
    "wood-rui": "0.2.0"
}

for pkg, ver in packages.items():
    if not check_version(pkg, ver):
        if pkg == "compas":
            uninstall_package(pkg)
        install_package(f"{pkg}=={ver}")

from compas_wood.binding_rhino import mesh_skeleton

if mesh:
    polylines = mesh_skeleton(mesh)

However, the results of the centerline detection are not very precise. Are there any settings in this library that can be adjusted to improve the accuracy of the output lines???

@Petras_Vestartas

I applied these settings, but they had no effect on the result:

#! python3
#venv: timber_connections
#r: compas==2.8.1, compas-wood==2.4.0, wood-nano==0.2.0, wood-rui==0.2.0

import subprocess
import sys

def check_version(pkg, target_version):
    try:
        import importlib.metadata as metadata
        version = metadata.version(pkg)
        return version == target_version
    except ImportError:
        return False
    except metadata.PackageNotFoundError:
        return False

def install_package(pkg):
    subprocess.check_call([sys.executable, "-m", "pip", "install", "--upgrade", pkg])

def uninstall_package(pkg):
    subprocess.call([sys.executable, "-m", "pip", "uninstall", "-y", pkg])

packages = {
    "compas": "2.8.1",
    "compas-wood": "2.4.0",
    "wood-nano": "0.2.0",
    "wood-rui": "0.2.0"
}

for pkg, ver in packages.items():
    if not check_version(pkg, ver):
        if pkg == "compas":
            uninstall_package(pkg)
        install_package(f"{pkg}=={ver}")

from compas_wood.binding_rhino import mesh_skeleton

if mesh:
    try:
        polylines = mesh_skeleton(
            mesh,
            voxel_size=0.03,
            smooth_iterations=12,
            tolerance=0.001,
            fit_nurbs=True
        )
    except TypeError:
        polylines = mesh_skeleton(mesh)

You need to have nicely remeshed closed triangular mesh.

Otherwise you can always check original c++ cgal implementation: