diff --git a/.github/workflows/build_wheels.yml b/.github/workflows/build_wheels.yml index 6ef0ce8..3704dec 100644 --- a/.github/workflows/build_wheels.yml +++ b/.github/workflows/build_wheels.yml @@ -14,30 +14,98 @@ on: - published env: - # Build `universal2` and `arm64` wheels on an Intel runner. - # Note that the `arm64` wheel and the `arm64` part of the `universal2` - # wheel cannot be tested in this configuration. CIBW_ARCHS_MACOS: "x86_64 arm64" + CIBW_ARCHS_WINDOWS: "AMD64" + CIBW_SKIP: "cp38-*" + CIBW_BUILD_VERBOSITY: 1 jobs: + # -------------------------------------------------------------------------- + # Build wheels for all OS × Python versions + # -------------------------------------------------------------------------- build_wheels: - name: Build wheels on ${{ matrix.os }} + name: Build wheels (${{ matrix.os }} / Python ${{ matrix.python-version }}) runs-on: ${{ matrix.os }} + strategy: + fail-fast: false matrix: - os: [ubuntu-latest, windows-latest, macos-latest] + include: + # ---------------- Ubuntu ---------------- + - os: ubuntu-latest + python-version: "3.9" + tag: cp39 + - os: ubuntu-latest + python-version: "3.10" + tag: cp310 + - os: ubuntu-latest + python-version: "3.11" + tag: cp311 + - os: ubuntu-latest + python-version: "3.12" + tag: cp312 + - os: ubuntu-latest + python-version: "3.13" + tag: cp313 + - os: ubuntu-latest + python-version: "3.14" + tag: cp314 + + # ---------------- Windows ---------------- + - os: windows-latest + python-version: "3.9" + tag: cp39 + - os: windows-latest + python-version: "3.10" + tag: cp310 + - os: windows-latest + python-version: "3.11" + tag: cp311 + - os: windows-latest + python-version: "3.12" + tag: cp312 + - os: windows-latest + python-version: "3.13" + tag: cp313 + - os: windows-latest + python-version: "3.14" + tag: cp314 + + # ---------------- macOS ---------------- + - os: macos-latest + python-version: "3.9" + tag: cp39 + - os: macos-latest + python-version: "3.10" + tag: cp310 + - os: macos-latest + python-version: "3.11" + tag: cp311 + - os: macos-latest + python-version: "3.12" + tag: cp312 + - os: macos-latest + python-version: "3.13" + tag: cp313 + - os: macos-latest + python-version: "3.14" + tag: cp314 + + env: + CIBW_BUILD: "${{ matrix.tag }}-*" steps: - uses: actions/checkout@v4 - name: Build wheels - uses: pypa/cibuildwheel@v3.2.1 + uses: pypa/cibuildwheel@v3.3.0 - uses: actions/upload-artifact@v4 with: - name: wheels-${{ matrix.os }} - path: ./wheelhouse/*.whl + name: wheels-${{ matrix.os }}-py${{ matrix.python-version }} + path: wheelhouse/*.whl + if-no-files-found: error build_sdist: name: Build source distribution diff --git a/.gitignore b/.gitignore index 56c814d..695379b 100644 --- a/.gitignore +++ b/.gitignore @@ -150,3 +150,4 @@ nose_debug test.ply .vscode *.ply +*.stl \ No newline at end of file diff --git a/README.md b/README.md index afe5c83..895fab1 100644 --- a/README.md +++ b/README.md @@ -2,17 +2,20 @@ [![CI](https://github.com/ifilot/pytessel/actions/workflows/build_wheels.yml/badge.svg)](https://github.com/ifilot/pytessel/actions/workflows/build_wheels.yml) [![PyPI](https://img.shields.io/pypi/v/pytessel?color=green)](https://pypi.org/project/pytessel/) +[![PyPI - Downloads](https://img.shields.io/pypi/dm/pypi)](https://pypi.org/project/pytessel/) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) ## Purpose -Python package for building isosurfaces from 3D scalar fields +PyTessel is a Python package for constructing isosurfaces from 3D scalar fields +using the marching cubes algorithm. It is designed for scientific visualization, +computational geometry, and mesh generation workflows. While PyTessel was +originally developed for rendering molecular orbitals, it is flexible enough to +tessellate arbitrary scalar fields. -## Installation +![isosurface](img/isosurface.png) -[![PyPI](https://img.shields.io/pypi/v/pytessel?color=green)](https://pypi.org/project/pytessel/) -[![PyPI - Downloads](https://img.shields.io/pypi/dm/pypi)](https://pypi.org/project/pytessel/) -![PyPI - Python Version](https://img.shields.io/pypi/pyversions/pytessel) +## Installation ``` pip install pytessel @@ -20,9 +23,13 @@ pip install pytessel ## Getting started -In the script below, the isosurface of a threedimensional Gaussian function is -constructed. The isosurface is written to `test.ply` and can, for example, -be opened using `ctmviewer` (Linux) or `3D viewer` (Windows). +The example below constructs an isosurface of a three-dimensional Gaussian function. +The resulting surface is written to `test.ply`, which can be viewed using tools +such as: + +* `ctmviewer` (Linux) +* `3D Viewer` (Windows, free via Microsoft Store) +* Blender, MeshLab, etc. ```python from pytessel import PyTessel @@ -31,21 +38,34 @@ import numpy as np def main(): pytessel = PyTessel() - # generate some data + # Generate a regular grid x = np.linspace(0, 10, 50) - # the grid is organized with z the slowest moving index and x the fastest moving index - grid = np.flipud(np.vstack(np.meshgrid(x, x, x, indexing='ij')).reshape(3,-1)).T - R = [5,5,5] - scalarfield = np.reshape(np.array([gaussian(r,R) for r in grid]), (len(x),len(x),len(x))) + # Grid ordering: + # z is the slowest-moving index, x is the fastest-moving index + grid = np.flipud( + np.vstack(np.meshgrid(x, x, x, indexing='ij')).reshape(3, -1) + ).T + + R = [5, 5, 5] + scalarfield = np.reshape( + np.array([gaussian(r, R) for r in grid]), + (len(x), len(x), len(x)) + ) + unitcell = np.diag(np.ones(3) * 10.0) - vertices, normals, indices = pytessel.marching_cubes(scalarfield.flatten(), scalarfield.shape, unitcell.flatten(), 0.1) + vertices, normals, indices = pytessel.marching_cubes( + scalarfield.flatten(), + scalarfield.shape, + unitcell.flatten(), + 0.1 + ) pytessel.write_ply('test.ply', vertices, normals, indices) def gaussian(r, R): - return np.exp(-(r-R).dot((r-R))) + return np.exp(-(r - R).dot(r - R)) if __name__ == '__main__': main() @@ -53,9 +73,17 @@ if __name__ == '__main__': ## Isosurface quality -In the script below, 6 different images are created of an icosahedral "metaball" using a grid -of `10x10x10`,`20x20x20`,`25x25x25`,`50x50x50`,`100x100x100`, and `200x200x200` points. The -resulting `.ply` files are rendered using [Blender](https://www.blender.org/). +The script below demonstrates how grid resolution affects surface quality. Six +isosurfaces of an icosahedral metaball are generated using grids of: + +* `10×10×10` +* `20×20×20` +* `25×25×25` +* `50×50×50` +* `100×100×100` +* `200×200×200` + +Each surface is exported as a .ply file and rendered using [Blender](https://www.blender.org/). ```python from pytessel import PyTessel @@ -67,50 +95,52 @@ def main(): """ pytessel = PyTessel() - for nrpoints in [10,20,25,50,100,200]: + for nrpoints in [10, 20, 25, 50, 100, 200]: sz = 3 x = np.linspace(-sz, sz, nrpoints) y = np.linspace(-sz, sz, nrpoints) z = np.linspace(-sz, sz, nrpoints) - xx, yy, zz, field = icosahedron_field(x,y,z) + xx, yy, zz, field = icosahedron_field(x, y, z) unitcell = np.diag(np.ones(3) * sz * 2) - pytessel = PyTessel() isovalue = 3.75 - vertices, normals, indices = pytessel.marching_cubes(field.flatten(), field.shape, unitcell.flatten(), isovalue) - pytessel.write_ply('icosahedron_%03i.ply' % nrpoints, vertices, normals, indices) - -def icosahedron_field(x,y,z): + vertices, normals, indices = pytessel.marching_cubes( + field.flatten(), + field.shape, + unitcell.flatten(), + isovalue + ) + + pytessel.write_ply( + f'icosahedron_{nrpoints:03d}.ply', + vertices, + normals, + indices + ) + +def icosahedron_field(x, y, z): """ - Produce a scalar field for the icosahedral metaballs + Produce a scalar field for icosahedral metaballs """ phi = (1 + np.sqrt(5)) / 2 vertices = [ - [0,1,phi], - [0,-1,-phi], - [0,1,-phi], - [0,-1,phi], - [1,phi,0], - [-1,-phi,0], - [1,-phi,0], - [-1,phi,0], - [phi,0,1], - [-phi,0,-1], - [phi,0,-1], - [-phi,0,1] + [0, 1, phi], [0, -1, -phi], [0, 1, -phi], [0, -1, phi], + [1, phi, 0], [-1, -phi, 0], [1, -phi, 0], [-1, phi, 0], + [phi, 0, 1], [-phi, 0, -1], [phi, 0, -1], [-phi, 0, 1] ] - xx,yy,zz = np.meshgrid(x,y,z) + xx, yy, zz = np.meshgrid(x, y, z) field = np.zeros_like(xx) + for v in vertices: - field += f(xx,yy,zz,v[0], v[1],v[2]) + field += metaball(xx, yy, zz, v[0], v[1], v[2]) - return xx,yy,zz,field + return xx, yy, zz, field -def f(x,y,z,X0,Y0,Z0): +def metaball(x, y, z, X0, Y0, Z0): """ Single metaball function """ @@ -120,10 +150,14 @@ if __name__ == '__main__': main() ``` -![Icosahedral metaballs](https://raw.githubusercontent.com/ifilot/pytessel/master/img/metaballs_3x2.png) +![Icosahedral metaballs](img/metaballs_3x2.png) -## Local compilation (Linux) +## Gallery -```bash -python3 setup.py build_ext --inplace bdist -``` +### Stanford Bunny + +![bunny](img/bunny.png) + +### Gyroid surface + +![gyroid](img/gyroid.png) \ No newline at end of file diff --git a/doc.Dockerfile b/doc.Dockerfile deleted file mode 100644 index 000052c..0000000 --- a/doc.Dockerfile +++ /dev/null @@ -1,2 +0,0 @@ -FROM nginx:latest -COPY html-docs/ /usr/share/nginx/html/ \ No newline at end of file diff --git a/environment.yml b/environment.yml deleted file mode 100644 index 92f994f..0000000 --- a/environment.yml +++ /dev/null @@ -1,6 +0,0 @@ -name: pytessel -dependencies: - - numpy - - conda-build - - conda-verify - - anaconda-client diff --git a/example/build_metaballs_icosahedron_rectangular_grid.py b/example/build_metaballs_icosahedron_rectangular_grid.py deleted file mode 100644 index 48e65b7..0000000 --- a/example/build_metaballs_icosahedron_rectangular_grid.py +++ /dev/null @@ -1,60 +0,0 @@ -from pytessel import PyTessel -import numpy as np - -def main(): - """ - Build 6 .ply files of increasing quality - """ - pytessel = PyTessel() - - sz = 3 - - x = np.linspace(-sz, sz, 30) - y = np.linspace(-sz, sz, 40) - z = np.linspace(-sz, sz, 50) - - xx, yy, zz, field = icosahedron_field(x,y,z) - print(field.shape) - - unitcell = np.diag(np.ones(3) * sz * 2) - pytessel = PyTessel() - isovalue = 3.75 - vertices, normals, indices = pytessel.marching_cubes(field.flatten(), field.shape, unitcell.flatten(), isovalue) - - pytessel.write_ply('icosahedron_metaballs.ply', vertices, normals, indices) - -def icosahedron_field(x,y,z): - """ - Produce a scalar field for the icosahedral metaballs - """ - phi = (1 + np.sqrt(5)) / 2 - vertices = [ - [0,1,phi], - [0,-1,-phi], - [0,1,-phi], - [0,-1,phi], - [1,phi,0], - [-1,-phi,0], - [1,-phi,0], - [-1,phi,0], - [phi,0,1], - [-phi,0,-1], - [phi,0,-1], - [-phi,0,1] - ] - - zz,yy,xx = np.meshgrid(z,y,x, indexing='ij') - field = np.zeros_like(xx) - for v in vertices: - field += f(xx,yy,zz,v[0], v[1],v[2]) - - return xx, yy, zz, field - -def f(x,y,z,X0,Y0,Z0): - """ - Single metaball function - """ - return 1 / ((x - X0)**2 + (y - Y0)**2 + (z - Z0)**2) - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/example/.gitignore b/examples/.gitignore similarity index 100% rename from example/.gitignore rename to examples/.gitignore diff --git a/examples/bunny.py b/examples/bunny.py new file mode 100644 index 0000000..d15d47f --- /dev/null +++ b/examples/bunny.py @@ -0,0 +1,57 @@ +import numpy as np +from scipy.ndimage import gaussian_filter +from pytessel import PyTessel +import os + +ROOT = os.path.dirname(__file__) + +def main(): + output_ply = os.path.join(ROOT, 'bunny.ply') + + # load data points + data = np.load(os.path.join(ROOT, 'bunny_pointcloud.npz')) + vertices = data["vertices"] + + print("Voxelizing...") + grid = voxelize(vertices, n=256) + + print("Smoothing...") + field = gaussian_filter(grid, sigma=2.0) + field /= field.max() # normalize + + print("Running marching cubes...") + t = PyTessel() + v, nrm, idx = t.marching_cubes( + field.ravel(), + dimensions=field.shape, + unitcell=(1, 0, 0, + 0, 1, 0, + 0, 0, 1), + isovalue=0.01, + ) + + print(f"Extracted {len(v)} vertices") + + print("Writing output...") + t.write_ply(output_ply, v, nrm, idx) + + print(f"Done! Output written to {output_ply}") + +def voxelize(points, n=256, padding=0.05): + mins = points.min(axis=0) + maxs = points.max(axis=0) + + size = maxs - mins + mins -= padding * size + maxs += padding * size + + grid = np.zeros((n, n, n), dtype=np.float32) + + idx = ((points - mins) / (maxs - mins) * (n - 1)).astype(np.int32) + idx = np.clip(idx, 0, n - 1) + + grid[idx[:, 0], idx[:, 1], idx[:, 2]] += 1.0 + return grid + +if __name__ == "__main__": + main() diff --git a/examples/bunny_pointcloud.npz b/examples/bunny_pointcloud.npz new file mode 100644 index 0000000..7a8f5e9 Binary files /dev/null and b/examples/bunny_pointcloud.npz differ diff --git a/examples/gyroid.py b/examples/gyroid.py new file mode 100644 index 0000000..5f8db46 --- /dev/null +++ b/examples/gyroid.py @@ -0,0 +1,42 @@ +import numpy as np +from pytessel import PyTessel +import os + +ROOT = os.path.dirname(__file__) + +def generate_gyroid( + n=192, + periods=2, +): + """ + Generate a gyroid scalar field on a cubic grid. + """ + x = np.linspace(0, 2 * np.pi * periods, n, endpoint=False) + X, Y, Z = np.meshgrid(x, x, x, indexing="ij") + + field = ( + np.sin(X) * np.cos(Y) + + np.sin(Y) * np.cos(Z) + + np.sin(Z) * np.cos(X) + ) + + return field.astype(np.float32) + + +if __name__ == "__main__": + t = PyTessel() + + n = 256 + grid = generate_gyroid(n=n, periods=2) + + vertices, normals, indices = t.marching_cubes( + grid.ravel(), + dimensions=(n, n, n), + unitcell=(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0), + isovalue=0.01, + ) + + t.write_ply(os.path.join(ROOT,"gyroid.ply"), vertices, normals, indices) + t.write_stl(os.path.join(ROOT, "gyroid.stl"), vertices, normals, indices) diff --git a/example/build_metaballs_icosahedron.py b/examples/metaballs_icosahedron.py similarity index 93% rename from example/build_metaballs_icosahedron.py rename to examples/metaballs_icosahedron.py index 9f3c95a..7a54cc1 100644 --- a/example/build_metaballs_icosahedron.py +++ b/examples/metaballs_icosahedron.py @@ -29,6 +29,7 @@ def main(): isovalue) pytessel.write_ply('icosahedron_%03i.ply' % nrpoints, vertices, normals, indices) + pytessel.write_stl('icosahedron_%03i.stl' % nrpoints, vertices, normals, indices) # perform an additional test wherein an isosurface is generated using a # field for which the grid dimensions differ per cartesian direction diff --git a/example/metaballs_scaling.py b/examples/metaballs_scaling.py similarity index 100% rename from example/metaballs_scaling.py rename to examples/metaballs_scaling.py diff --git a/img/bunny.png b/img/bunny.png new file mode 100644 index 0000000..30f2a5d Binary files /dev/null and b/img/bunny.png differ diff --git a/img/gyroid.png b/img/gyroid.png new file mode 100644 index 0000000..33b0691 Binary files /dev/null and b/img/gyroid.png differ diff --git a/img/isosurface.png b/img/isosurface.png new file mode 100644 index 0000000..9fc920b Binary files /dev/null and b/img/isosurface.png differ diff --git a/img/metaballs_3x2.png b/img/metaballs_3x2.png index 86ee936..44d2695 100644 Binary files a/img/metaballs_3x2.png and b/img/metaballs_3x2.png differ diff --git a/meson.build b/meson.build new file mode 100644 index 0000000..580ecc0 --- /dev/null +++ b/meson.build @@ -0,0 +1,103 @@ +# --------------------------------------------------------------------------- +# Top-level project definition +# --------------------------------------------------------------------------- +# +# This project builds a Python extension module backed by C++ and Cython. +# We explicitly require a modern Meson version to ensure stable behavior +# with meson-python and PEP 517 builds. +# +project( + 'pytessel', + ['cpp', 'cython'], + meson_version: '>=1.0.0', + default_options: [ + 'buildtype=release', # Optimized builds by default + 'cpp_std=c++17', # Required by the C++ sources + ] +) + +# --------------------------------------------------------------------------- +# Python and compiler discovery +# --------------------------------------------------------------------------- +# +# We ask Meson for the *non-pure* Python installation because this project +# builds a compiled extension (.pyd / .so), not a pure-Python wheel. +# +# meson-python ensures that this Python installation corresponds to the +# wheel being built inside the PEP 517 environment. +# +py = import('python').find_installation(pure: false) + +# Fetch the C++ compiler object. We do not force a specific compiler here; +# the actual toolchain selection is handled externally (CI configuration). +# +# On Windows CI we intentionally rely on MinGW for robustness. +# +cpp = meson.get_compiler('cpp') + + +# --------------------------------------------------------------------------- +# OpenMP support +# --------------------------------------------------------------------------- +openmp_dep = [] +openmp = dependency('openmp', required: false) +if openmp.found() + openmp_dep = [openmp] +endif + +# --------------------------------------------------------------------------- +# Include directories +# --------------------------------------------------------------------------- +# +# Project headers live inside the Python package directory itself. +# This allows: +# - a clean package layout +# - easy access to headers from both C++ and Cython sources +# +inc = include_directories('pytessel') + + +# --------------------------------------------------------------------------- +# Compiled Python extension module +# --------------------------------------------------------------------------- +# +# This builds the core extension module (pytessel_core) which is imported from +# Python as: +# +# from pytessel import pytessel_core +# +# Notes: +# - Sources include both C++ and a Cython .pyx file +# - cython_language=cpp ensures Cython generates C++ code +# - subdir='pytessel' ensures the extension is installed inside the package +# - install=true is required for wheel builds +# +# Toolchain notes (important for CI): +# - On Windows, this is built with MinGW in CI and repaired with delvewheel +# - We intentionally avoid MSVC-specific logic here to keep meson.build +# platform-agnostic and CI-stable +# +py.extension_module( + 'pytessel_core', + [ + 'pytessel/pytessel_core.pyx', + 'pytessel/isosurface_mesh.cpp', + 'pytessel/isosurface.cpp', + 'pytessel/scalar_field.cpp', + ], + subdir: 'pytessel', + include_directories: inc, + override_options: ['cython_language=cpp'], + dependencies: openmp_dep, + install: true, +) + +# --------------------------------------------------------------------------- +# Pure-Python source files +# --------------------------------------------------------------------------- +py.install_sources( + 'pytessel/__init__.py', + 'pytessel/_version.py', + + subdir: 'pytessel', +) \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 7834a23..3f7458f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,15 +1,131 @@ +# --------------------------------------------------------------------------- +# Build system configuration +# --------------------------------------------------------------------------- +# +# We use meson-python (PEP 517 backend) to build the C/C++ extension. +# In principle this works well across platforms, but Windows deserves +# special attention (see cibuildwheel.windows below). +# +# NOTE: +# - meson-python runs builds in a fully isolated PEP 517 environment. +# - Environment setup done in GitHub Actions (PATH, vcvarsall, etc.) +# does NOT reliably propagate into the actual Meson invocation. +# +# This isolation is the root cause of most Windows-specific complexity. +# [build-system] requires = [ - "setuptools>=42", - "wheel", + "meson-python", "Cython", "numpy", + + # Required only for Windows wheels when building with MinGW. + # delvewheel is used to bundle MinGW runtime DLLs into the wheel. + # + # IMPORTANT: + # We intentionally do NOT rely on MSVC here, because: + # - MSVC + Meson + meson-python + cibuildwheel is currently brittle + # in CI due to PATH/link.exe collisions with Git for Windows. + # - meson-python does not support platform-scoped arguments, making + # it impossible to reliably force MSVC via Meson native files + # without breaking Linux/macOS builds. + # + "delvewheel" +] +build-backend = "mesonpy" + + +# --------------------------------------------------------------------------- +# Project metadata +# --------------------------------------------------------------------------- +[project] +name = "pytessel" +version = "1.3.0" +description = "Python package for building isosurfaces from 3D scalar fields" +readme = "README.md" +license = { text = "GPL-3.0-only" } +authors = [ + { name = "Ivo Filot", email = "ivo@ivofilot.nl" } ] -build-backend = "setuptools.build_meta" +requires-python = ">=3.9" + +# Runtime Python dependencies (independent of build backend) +dependencies = ["numpy"] + +[project.urls] +Homepage = "https://github.com/ifilot/pyqint" + +# --------------------------------------------------------------------------- +# cibuildwheel configuration +# --------------------------------------------------------------------------- +# +# cibuildwheel is used to build wheels for Linux, macOS, and Windows. +# +# IMPORTANT DESIGN DECISION (Windows): +# +# We intentionally build Windows wheels using MinGW rather than MSVC. +# +# Why? +# ---- +# Using MSVC with Meson in CI looks appealing, but in practice it is +# extremely fragile due to: +# +# - MSVC tool names (link.exe, lib.exe) colliding with Git for Windows +# - PATH ordering being rebuilt inside PEP 517 isolation +# - meson-python lacking platform-scoped configuration hooks +# - Meson correctly refusing mixed toolchains +# +# After extensive attempts, the MinGW + delvewheel approach proved to be: +# - deterministic +# - CI-stable +# - well-supported by cibuildwheel +# - transparent to end users +# +# The cost is a small number of bundled runtime DLLs, which delvewheel +# handles automatically. +# [tool.cibuildwheel] -test-requires = "pytest numpy cython" + +# pytest is used to validate the built wheel after installation +test-requires = "pytest" test-command = "pytest {project}/tests" -# skip PyPy wheels -skip = ["pp*", "cp36-*"] +# Skip configurations we do not support or do not need +skip = [ + "pp*", # PyPy not supported + "*-win32", # 32-bit Windows not supported + "*-manylinux_i686", # 32-bit Linux not supported + "cp*-musllinux_*", # musllinux not supported + "cp38-*", # Python < 3.9 not supported +] + + +# --------------------------------------------------------------------------- +# Windows-specific wheel repair +# --------------------------------------------------------------------------- +# +# MinGW-built extensions depend on runtime DLLs such as: +# - libstdc++-6.dll +# - libgcc_s_seh-1.dll +# - libwinpthread-1.dll +# +# These DLLs are NOT present on a default Windows system. +# +# delvewheel: +# - scans the built .pyd file +# - finds MinGW runtime dependencies +# - bundles them into the wheel +# - patches the loader so imports work out-of-the-box +# +# This step is REQUIRED for functional Windows wheels when using MinGW. +# +[tool.cibuildwheel.windows] + +# Install delvewheel into the isolated cibuildwheel environment +before-build = [ + "pip install delvewheel", +] + +# Repair the wheel after build by bundling MinGW runtime DLLs +repair-wheel-command = "delvewheel repair -w {dest_dir} {wheel}" \ No newline at end of file diff --git a/pytessel/__init__.py b/pytessel/__init__.py index b089e54..edea36e 100644 --- a/pytessel/__init__.py +++ b/pytessel/__init__.py @@ -1,3 +1,3 @@ -from .pytessel import PyTessel +from .pytessel_core import PyTessel from ._version import __version__ diff --git a/pytessel/_version.py b/pytessel/_version.py old mode 100755 new mode 100644 index 7f647d0..d46ba33 --- a/pytessel/_version.py +++ b/pytessel/_version.py @@ -1 +1,6 @@ -__version__ = "1.2.7" +from importlib.metadata import version, PackageNotFoundError + +try: + __version__ = version("pytessel") +except PackageNotFoundError: + __version__ = "0.0.0" # optional fallback \ No newline at end of file diff --git a/pytessel/pytessel.pyx b/pytessel/pytessel.pyx deleted file mode 100644 index 9cc9dbb..0000000 --- a/pytessel/pytessel.pyx +++ /dev/null @@ -1,135 +0,0 @@ -# distutils: language = c++ -# cython: c_string_type=unicode, c_string_encoding=utf8 - -from .pytessel cimport ScalarField, IsoSurface -from libcpp.string cimport string -from libcpp.memory cimport shared_ptr,make_shared -import numpy as np -import sys -import cython -import numpy.typing as npt - -cdef class PyTessel: - - def __cinit__(self): - pass - - @cython.embedsignature(True) - def marching_cubes(self, vector[float] grid, vector[size_t] dimensions, vector[float] unitcell, float isovalue) -> \ - (npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64]): - """ - Perform marching cubes algorithm to generate isosurface - - Parameters - ---------- - grid : Iterable of floats - Scalar field as a flattened array - dimensions : Iterable of ints - Dimensions of the scalar field grid (nx, ny, nz) - unitcell : Iterable of floats - Unitcell matrix (flattened) - isovalue : float - Isovalue of the isosurface - - Returns - ------- - vertices : (Nx3) numpy array of floats - Triangle vertices - normals : (Nx3) numpy array of floats - Triangle normals (at the vertices) - indices : numpy array of ints - Triangle indices - - Notes - ----- - * The scalar field needs to be encoded such that the z-coordinate is the slowest moving - index and the x-coordinate the fastest moving index. - * The dimensions of the scalar field are encoded in the order :code:`(nx, ny, nz)`. - * You can use :code:`reversed(grid.shape)` to pass the dimensions. - * The output of the :code:`marching_cubes` function are three arrays, corresponding to the - vertices, the normals, and the indices of the isosurface. The isosurface corresponds to - a number of linked triangles (polygons) whose vertices are stored in the :code:`vertices` - array. For each vertex, the normal vector is encoded in the :code:`normals` array. Finally, - the triangles are stored as a triplet of indices in the :code:`indices` array. These indices - refer to the position in the :code:`vertices` and :code:`normals` array. Because multiple - triangles can use the same vertices, this is an efficient way to store the isosurface. - * One rarely needs to perform any operations on the :code:`vertices`, :code:`normals` and - :code:`indices` arrays. Typically, these arrays are constructed and immediately relayed - to the :code:`write_ply` function to store them as a file which can be used in another - program. - """ - cdef shared_ptr[ScalarField] scalarfield - cdef shared_ptr[IsoSurface] isosurface - cdef shared_ptr[IsoSurfaceMesh] isosurface_mesh - - # build scalar field - scalarfield = make_shared[ScalarField](grid, dimensions, unitcell) - - # construct isosurface - isosurface = make_shared[IsoSurface](scalarfield) - isosurface.get().marching_cubes(isovalue) - - # extract isosurface mesh - isosurface_mesh = make_shared[IsoSurfaceMesh](scalarfield, isosurface) - isosurface_mesh.get().construct_mesh(False) - - # extract data - vertices = np.array(isosurface_mesh.get().get_vertices(), dtype=np.float32).reshape(-1,3) - normals = np.array(isosurface_mesh.get().get_normals(), dtype=np.float32).reshape(-1,3) - indices = np.array(isosurface_mesh.get().get_indices(), dtype=np.uint32) - - return vertices, normals, indices - - @cython.embedsignature(True) - def write_ply(self, filename:str, vertices:npt.NDArray[np.float64], normals:npt.NDArray[np.float64], indices:npt.NDArray[np.uint32]) -> None: - """Stores isosurface as .ply file - - Parameters - ---------- - filename : str - Path to file - vertices : Iterable of floats - Array of vertices - normals : Iterable of floats - Array of normals - indices : Iterable of ints - Array of triangle indices - - Returns - ------- - None - - Notes - ----- - * This function is designed to be used with :code:`marching_cubes`. One can directly relay the output - of :code:`marching_cube` as the input for :code:`write_ply`. - """ - f = open(filename, 'wb') - - f.write(b"ply\n") - if sys.byteorder == 'little': - f.write(b"format binary_little_endian 1.0\n") - else: - f.write(b"format binary_big_endian 1.0\n") - - f.write(b"comment test\n") - f.write(b"element vertex %i\n" % vertices.shape[0]) - f.write(b"property float x\n") - f.write(b"property float y\n") - f.write(b"property float z\n") - f.write(b"property float nx\n") - f.write(b"property float ny\n") - f.write(b"property float nz\n") - f.write(b"element face %i\n" % int(len(indices)/3)) - f.write(b"property list uchar uint vertex_indices\n") - f.write(b"end_header\n") - - for i in range(0, len(vertices)): - f.write(vertices[i].tobytes()) - f.write(normals[i].tobytes()) - - for i in range(0, int(len(indices)/3)): - f.write(b'\x03') - f.write(indices[i*3:(i+1)*3].tobytes()) - - f.close() diff --git a/pytessel/pytessel.pxd b/pytessel/pytessel_core.pxd similarity index 80% rename from pytessel/pytessel.pxd rename to pytessel/pytessel_core.pxd index 064da3d..00f0185 100644 --- a/pytessel/pytessel.pxd +++ b/pytessel/pytessel_core.pxd @@ -1,21 +1,12 @@ -# distutils: language = c++ +# cython: language_level=3 +# cython: language=c++ +# cython: module_name=pytessel_core +# cython: c_string_type=unicode, c_string_encoding=utf8 from libcpp.vector cimport vector from libcpp.string cimport string from libcpp.memory cimport shared_ptr -# IsoSurface -cdef extern from "isosurface.cpp": - pass - -# Isosurface Mesh -cdef extern from "isosurface_mesh.cpp": - pass - -# Scalar Field -cdef extern from "scalar_field.cpp": - pass - # Scalar Field class cdef extern from "scalar_field.h": cdef cppclass ScalarField: diff --git a/pytessel/pytessel_core.pyx b/pytessel/pytessel_core.pyx new file mode 100644 index 0000000..01a7818 --- /dev/null +++ b/pytessel/pytessel_core.pyx @@ -0,0 +1,215 @@ +# cython: language_level=3 +# cython: language=c++ +# cython: module_name=pytessel_core +# cython: c_string_type=unicode, c_string_encoding=utf8 + +from .pytessel_core cimport ScalarField, IsoSurface +from libcpp.string cimport string +from libcpp.memory cimport shared_ptr,make_shared +import numpy as np +import sys +import cython +import numpy.typing as npt + +cdef class PyTessel: + + def __cinit__(self): + pass + + @cython.embedsignature(True) + def marching_cubes( + self, + vector[float] grid, + vector[size_t] dimensions, + vector[float] unitcell, + float isovalue + ) -> tuple[ + npt.NDArray[np.float64], + npt.NDArray[np.float64], + npt.NDArray[np.float64] + ]: + """ + Perform marching cubes algorithm to generate isosurface + + Parameters + ---------- + grid : Iterable of floats + Scalar field as a flattened array + dimensions : Iterable of ints + Dimensions of the scalar field grid (nx, ny, nz) + unitcell : Iterable of floats + Unitcell matrix (flattened) + isovalue : float + Isovalue of the isosurface + + Returns + ------- + vertices : (Nx3) numpy array of floats + Triangle vertices + normals : (Nx3) numpy array of floats + Triangle normals (at the vertices) + indices : numpy array of ints + Triangle indices + + Notes + ----- + * The scalar field needs to be encoded such that the z-coordinate is the slowest moving + index and the x-coordinate the fastest moving index. + * The dimensions of the scalar field are encoded in the order :code:`(nx, ny, nz)`. + * You can use :code:`reversed(grid.shape)` to pass the dimensions. + * The output of the :code:`marching_cubes` function are three arrays, corresponding to the + vertices, the normals, and the indices of the isosurface. The isosurface corresponds to + a number of linked triangles (polygons) whose vertices are stored in the :code:`vertices` + array. For each vertex, the normal vector is encoded in the :code:`normals` array. Finally, + the triangles are stored as a triplet of indices in the :code:`indices` array. These indices + refer to the position in the :code:`vertices` and :code:`normals` array. Because multiple + triangles can use the same vertices, this is an efficient way to store the isosurface. + * One rarely needs to perform any operations on the :code:`vertices`, :code:`normals` and + :code:`indices` arrays. Typically, these arrays are constructed and immediately relayed + to the :code:`write_ply` function to store them as a file which can be used in another + program. + """ + cdef shared_ptr[ScalarField] scalarfield + cdef shared_ptr[IsoSurface] isosurface + cdef shared_ptr[IsoSurfaceMesh] isosurface_mesh + + # build scalar field + scalarfield = make_shared[ScalarField](grid, dimensions, unitcell) + + # construct isosurface + isosurface = make_shared[IsoSurface](scalarfield) + isosurface.get().marching_cubes(isovalue) + + # extract isosurface mesh + isosurface_mesh = make_shared[IsoSurfaceMesh](scalarfield, isosurface) + isosurface_mesh.get().construct_mesh(False) + + # extract data + vertices = np.array(isosurface_mesh.get().get_vertices(), dtype=np.float32).reshape(-1,3) + normals = np.array(isosurface_mesh.get().get_normals(), dtype=np.float32).reshape(-1,3) + indices = np.array(isosurface_mesh.get().get_indices(), dtype=np.uint32) + + return vertices, normals, indices + + def write_ply(self, + filename: str, + vertices: npt.NDArray[np.float64], + normals: npt.NDArray[np.float64], + indices: npt.NDArray[np.uint32], + ) -> None: + """ + Write a binary PLY file with vertices, normals, and triangular faces. + """ + + if vertices.shape != normals.shape: + raise ValueError("vertices and normals must have the same shape") + + if vertices.shape[1] != 3: + raise ValueError("vertices must be of shape (N, 3)") + + if indices.ndim != 1 or len(indices) % 3 != 0: + raise ValueError("indices must be a flat array of length multiple of 3") + + n_vertices = vertices.shape[0] + n_faces = len(indices) // 3 + + endian = "binary_little_endian" if sys.byteorder == "little" else "binary_big_endian" + + header = ( + "ply\n" + f"format {endian} 1.0\n" + "comment generated by write_ply\n" + f"element vertex {n_vertices}\n" + "property float x\n" + "property float y\n" + "property float z\n" + "property float nx\n" + "property float ny\n" + "property float nz\n" + f"element face {n_faces}\n" + "property list uchar uint vertex_indices\n" + "end_header\n" + ).encode("ascii") + + # Interleave vertices and normals: (N, 6) + vertex_data = np.hstack((vertices, normals)).astype(np.float32, copy=False) + + # Face data: [3, i0, i1, i2] + face_data = np.empty((n_faces, 4), dtype=np.uint32) + face_data[:, 0] = 3 + face_data[:, 1:] = indices.reshape(-1, 3) + + # The first element must be uint8, the rest uint32 + face_bytes = np.empty( + n_faces, + dtype=[("n", "u1"), ("i", "u4", (3,))] + ) + face_bytes["n"] = 3 + face_bytes["i"] = indices.reshape(-1, 3) + + with open(filename, "wb") as f: + f.write(header) + f.write(vertex_data.tobytes()) + f.write(face_bytes.tobytes()) + + def write_stl(self, + filename: str, + vertices: npt.NDArray[np.float64], + normals: npt.NDArray[np.float64], + indices: npt.NDArray[np.uint32], + ) -> None: + """ + Write a binary STL file. + + Parameters + ---------- + vertices : (N, 3) array of vertex positions + indices : (M,) flat array, length multiple of 3 + normals : (N, 3) vertex normals + """ + + if vertices.shape != normals.shape: + raise ValueError("vertices and normals must have the same shape") + + if vertices.shape[1] != 3: + raise ValueError("vertices must be of shape (N, 3)") + + if indices.ndim != 1 or len(indices) % 3 != 0: + raise ValueError("indices must be a flat array with length multiple of 3") + + # Build triangle vertices: (T, 3, 3) + triangles = vertices[indices].reshape(-1, 3, 3) + + # Compute face normals from vertex normals + tri_normals = normals[indices].reshape(-1, 3, 3) + face_normals = tri_normals.sum(axis=1) + + # Normalize + norm = np.linalg.norm(face_normals, axis=1) + norm[norm == 0.0] = 1.0 + face_normals /= norm[:, None] + + n_triangles = triangles.shape[0] + + header = b"Generated by write_stl".ljust(80, b"\0") + + with open(filename, "wb") as f: + f.write(header) + f.write(np.uint32(n_triangles).tobytes()) + + dtype = np.dtype([ + ("normal", "f4", (3,)), + ("v1", "f4", (3,)), + ("v2", "f4", (3,)), + ("v3", "f4", (3,)), + ("attr", "u2"), + ]) + + data = np.empty(n_triangles, dtype=dtype) + data["normal"] = face_normals.astype(np.float32) + data["v1"] = triangles[:, 0].astype(np.float32) + data["v2"] = triangles[:, 1].astype(np.float32) + data["v3"] = triangles[:, 2].astype(np.float32) + data["attr"] = 0 + + f.write(data.tobytes()) diff --git a/pytessel/vec3.h b/pytessel/vec3.h index bb700f3..6319be6 100644 --- a/pytessel/vec3.h +++ b/pytessel/vec3.h @@ -1,6 +1,8 @@ #ifndef _VEC3_H #define _VEC3_H +#include + typedef float mat33[3][3]; /** diff --git a/setup.py b/setup.py deleted file mode 100755 index 36485cf..0000000 --- a/setup.py +++ /dev/null @@ -1,120 +0,0 @@ -from setuptools import Extension, setup -from Cython.Build import cythonize -import os -import sys -import re - -PKG = "pytessel" -VERSIONFILE = os.path.join(os.path.dirname(__file__), PKG, "_version.py") -verstr = "unknown" -try: - verstrline = open(VERSIONFILE, "rt").read() -except EnvironmentError: - pass # Okay, there is no version file. -else: - VSRE = r"^__version__ = ['\"]([^'\"]*)['\"]" - mo = re.search(VSRE, verstrline, re.M) - if mo: - verstr = mo.group(1) - else: - print(r"Unable to find version in %s" % (VERSIONFILE,)) - raise RuntimeError(r"If %s.py exists, it is required to be well-formed" % (VERSIONFILE,)) - -def find_windows_versions(): - """ - Autofind the msvc and winkit versions; this is mainly used for local development / compilation - """ - root = os.path.join('C:', os.sep,'Program Files', 'Microsoft Visual Studio', '2022', 'Community', 'VC', 'Tools', 'MSVC') - - # for Gitlab actions, the above folder does not exist and this is communicated - # back by providing None as the result - if not os.path.exists(root): - return None, None - - for file in os.listdir(root): - if os.path.isdir(os.path.join(root, file)): - msvcver = file - - root = os.path.join('C:', os.sep,'Program Files (x86)', 'Windows Kits', '10', 'Include') - for file in os.listdir(root): - if os.path.isdir(os.path.join(root, file)): - winkitver = file - - return msvcver, winkitver - -# specify paths on Windows to find compiler and libraries -if os.name == 'nt': - msvc_ver, winkit_ver = find_windows_versions() - - if msvc_ver and winkit_ver: - # only proceed with setting the paths for local development, i.e. when the - # msvc_ver and winkit_ver variables are *not* None - os.environ['PATH'] += r";C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\%s\bin\Hostx64\x64" % msvc_ver - os.environ['PATH'] += r";C:\Program Files (x86)\Windows Kits\10\bin\%s\x64" % winkit_ver - - # set path to include folders - os.environ['INCLUDE'] += r";C:\Program Files (x86)\Windows Kits\10\Include\%s\ucrt" % winkit_ver - os.environ['INCLUDE'] += r";C:\Program Files (x86)\Windows Kits\10\Include\%s\shared" % winkit_ver - os.environ['INCLUDE'] += r";C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\%s\include" % msvc_ver - - # some references to libraries - os.environ['LIB'] += r";C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\%s\lib\x64" % msvc_ver - os.environ['LIB'] += r";C:\Program Files (x86)\Windows Kits\10\Lib\%s\um\x64" % winkit_ver - os.environ['LIB'] += r";C:\Program Files (x86)\Windows Kits\10\Lib\%s\ucrt\x64" % winkit_ver - else: - # re-order paths to ensure that the MSVC toolchain is in front; this needs to be done - # because the Git bin folder precedes the MSVC bin folder, resulting in the wrong link.exe - # executable to be used in the linking step - paths = os.environ['PATH'].split(";") - newpaths = [] - for path in paths: - if "Microsoft Visual Studio" in path: - newpaths = [path] + newpaths - else: - newpaths.append(path) - os.environ['PATH'] = ";".join(newpaths) - -if os.name == 'posix' and sys.platform != 'darwin': - extra_compile_args = ["-Wno-date-time", "-fopenmp", "-fPIC"] - extra_link_args = ["-fopenmp"] -elif os.name == 'nt': - extra_compile_args = ["/wd4244"] - extra_link_args = [] -elif sys.platform == 'darwin': - extra_compile_args = ["-Wno-date-time", "-fPIC", "-std=c++14"] - extra_link_args = [] - -ext_modules = [ - Extension( - "pytessel.pytessel", - ["pytessel/pytessel.pyx"], - extra_compile_args=extra_compile_args, # overrule some arguments - extra_link_args=extra_link_args - ), -] - -with open("README.md", "r", encoding="utf-8") as fh: - long_description = fh.read() - -setup( - name=PKG, - version=verstr, - author="Ivo Filot", - author_email="ivo@ivofilot.nl", - description="Python package for building isosurfaces from 3D scalar fields", - long_description=long_description, - long_description_content_type="text/markdown", - url="https://github.com/ifilot/pytessel", - ext_modules=cythonize(ext_modules[0], - language_level = "3", - build_dir="build"), - packages=[PKG], - include_package_data=True, - classifiers=[ - "Programming Language :: Python :: 3", - "License :: OSI Approved :: MIT License", - "Operating System :: POSIX", - ], - python_requires='>=3.5', - install_requires=['numpy'], -) diff --git a/testversion.py b/testversion.py deleted file mode 100644 index a6f74ca..0000000 --- a/testversion.py +++ /dev/null @@ -1,72 +0,0 @@ -# -*- coding: utf-8 -*- - -import os -import re - -ROOT = os.path.dirname(__file__) - -def main(): - version_versionpy = get_version_versionpy() - version_metayaml = get_version_metayaml() - - print('Version strings found:') - print(version_versionpy) - print(version_metayaml) - - try: - for i in range(0,3): - assert version_versionpy[i] == version_metayaml[i] - except Exception as e: - print(e) - raise Exception('Invalid version strings encountered') - -def get_version_projecttoml(): - """ - Extract the version string from the pyproject.toml file - """ - pattern = re.compile(r'^version\s*=\s*"(\d+\.\d+.\d+)"\s*$') - - f = open(os.path.join(ROOT, 'pyproject.toml')) - lines = f.readlines() - for line in lines: - match = re.match(pattern, line) - if match: - version = match.groups(1)[0] - return [int(i) for i in version.split('.')] - - return None - -def get_version_versionpy(): - """ - Extract the version string from the _version.py file - """ - pattern = re.compile(r'^__version__\s*=\s*[\'"](\d+\.\d+.\d+)[\'"]\s*$') - - f = open(os.path.join(ROOT, 'pytessel', '_version.py')) - lines = f.readlines() - for line in lines: - match = re.match(pattern, line) - if match: - version = match.groups(1)[0] - return [int(i) for i in version.split('.')] - - return None - -def get_version_metayaml(): - """ - Extract the version string from the meta.yaml file - """ - pattern = re.compile(r'^\s*version\s*:\s*"(\d+\.\d+.\d+)"\s*$') - - f = open(os.path.join(ROOT, 'meta.yaml')) - lines = f.readlines() - for line in lines: - match = re.match(pattern, line) - if match: - version = match.groups(1)[0] - return [int(i) for i in version.split('.')] - - return None - -if __name__ == '__main__': - main()