diff --git a/.github/workflows/build-and-deploy.yml b/.github/workflows/build-and-deploy.yml index 269dabd..0348f24 100644 --- a/.github/workflows/build-and-deploy.yml +++ b/.github/workflows/build-and-deploy.yml @@ -24,7 +24,7 @@ jobs: - uses: actions/checkout@v4 - name: Build wheels - uses: pypa/cibuildwheel@v2.19.2 + uses: pypa/cibuildwheel@v2.22.0 - uses: actions/upload-artifact@v4 with: @@ -39,7 +39,7 @@ jobs: - uses: actions/setup-python@v5 with: - python-version: '3.12' + python-version: '3.13' - name: Build sdist run: pipx run build --sdist @@ -72,7 +72,7 @@ jobs: - name: List artifacts run: ls -R - - uses: pypa/gh-action-pypi-publish@v1.9.0 + - uses: pypa/gh-action-pypi-publish@v1.12.2 with: user: __token__ password: ${{ secrets.PYPI_TOKEN }} diff --git a/.gitignore b/.gitignore index ea381a1..3e82ccc 100644 --- a/.gitignore +++ b/.gitignore @@ -13,4 +13,5 @@ CMakeCache.txt CMakeFiles/ Makefile cmake_install.cmake -compile_commands.json \ No newline at end of file +compile_commands.json +.cache/ \ No newline at end of file diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 1878fb7..752950d 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -5,7 +5,7 @@ ci: repos: - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.5.2 + rev: v0.8.0 hooks: - id: ruff args: [--fix, --exit-non-zero-on-fix] @@ -32,14 +32,14 @@ repos: files: ^src/stl_reader/.*\.py - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.6.0 + rev: v5.0.0 hooks: - id: check-merge-conflict - id: debug-statements - id: trailing-whitespace - repo: https://github.com/pre-commit/mirrors-mypy - rev: v1.10.1 + rev: v1.13.0 hooks: - id: mypy additional_dependencies: [numpy>=2.0, npt-promote==0.1] @@ -54,12 +54,12 @@ repos: args: [--autofix, --indent, '2'] - repo: https://github.com/pre-commit/mirrors-clang-format - rev: v17.0.6 + rev: v19.1.4 hooks: - id: clang-format - repo: https://github.com/python-jsonschema/check-jsonschema - rev: 0.29.0 + rev: 0.29.4 hooks: - id: check-github-workflows diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 1bad15a..06db584 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -9,7 +9,7 @@ extensions. If using emacs and helm, generate the project configuration files using `-DCMAKE_EXPORT_COMPILE_COMMANDS=ON`. Here's a sample configuration for C++11 on Linux: ``` -pip install nanobind +pip install nanobind -U export NANOBIND_INCLUDE=$(python -c "import nanobind, os; print(os.path.join(os.path.dirname(nanobind.__file__), 'cmake'))") cmake -Dnanobind_DIR=$NANOBIND_INCLUDE -DCMAKE_EXPORT_COMPILE_COMMANDS=ON -DCMAKE_CXX_STANDARD_INCLUDE_DIRECTORIES="/usr/include/c++/11;/usr/include/x86_64-linux-gnu/c++/11/;/usr/lib/gcc/x86_64-linux-gnu/11/include/" ``` diff --git a/README.rst b/README.rst index a9e9259..1ae06cc 100644 --- a/README.rst +++ b/README.rst @@ -10,9 +10,11 @@ .. |MIT| image:: https://img.shields.io/badge/License-MIT-yellow.svg :target: https://opensource.org/licenses/MIT -``stl-reader`` is a Python library for raipidly reading binary STL -files. It wraps a Cython interface to the fast STL library provided by -`libstl `_. Thanks @aki5! +``stl-reader`` is a Python library for rapidly reading binary and ASCII +STL files. It wraps a `nanobind +`_ interface to the fast STL +library provided by `libstl `_. Thanks +@aki5! The main advantage of ``stl-reader`` over other STL reading libraries is its performance. It is particularly well-suited for large files, mainly @@ -167,6 +169,40 @@ Here's an additional benchmark comparing VTK with ``stl-reader``: .. image:: https://github.com/pyvista/stl-reader/raw/main/bench1.png +Read in ASCII Meshes +==================== + +The `stl-reader` also supports ASCII files and is around 2.4 times +faster than VTK at reading ASCII files. + +.. code:: python + + import time + import stl_reader + import pyvista as pv + import numpy as np + + # create and save a ASCII file + n = 1000 + mesh = pv.Plane(i_resolution=n, j_resolution=n).triangulate() + mesh.save("/tmp/tmp-ascii.stl", binary=False) + + # stl reader + tstart = time.time() + mesh = stl_reader.read_as_mesh("/tmp/tmp-ascii.stl") + print("stl-reader ", time.time() - tstart) + + tstart = time.time() + pv_mesh = pv.read("/tmp/tmp-ascii.stl") + print("pyvista reader", time.time() - tstart) + + # verify meshes are identical + assert np.allclose(mesh.points, pv_mesh.points) + + # approximate time to read in the 1M point file: + # stl-reader 0.80303955078125 + # pyvista reader 1.916085958480835 + ***************************** License and Acknowledgments ***************************** diff --git a/pyproject.toml b/pyproject.toml index 2a607bc..2c31eeb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -11,11 +11,11 @@ classifiers = [ "Operating System :: Microsoft :: Windows", "Operating System :: POSIX", "Intended Audience :: Science/Research", - "Programming Language :: Python :: 3.8", "Programming Language :: Python :: 3.9", "Programming Language :: Python :: 3.10", "Programming Language :: Python :: 3.11", - "Programming Language :: Python :: 3.12" + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13" ] dependencies = [ "numpy" @@ -23,12 +23,12 @@ dependencies = [ description = "Read in STL files" name = "stl-reader" readme = "README.rst" -requires-python = ">=3.8" -version = "0.2.1" +requires-python = ">=3.9" +version = "0.2.2" [tool.cibuildwheel] archs = ["auto64"] # 64-bit only -skip = "cp36-* cp37-* pp* *musllinux*" # disable PyPy and musl-based wheels +skip = "cp37-* pp* *musllinux*" # disable PyPy and musl-based wheels test-command = "pytest {project}/tests" test-requires = "pytest pyvista" diff --git a/src/stl_reader/_stlfile_wrapper.pyx b/src/stl_reader/_stlfile_wrapper.pyx deleted file mode 100644 index 0ae2bfc..0000000 --- a/src/stl_reader/_stlfile_wrapper.pyx +++ /dev/null @@ -1,92 +0,0 @@ -# cython: language_level=3 -# cython: boundscheck=False -# cython: wraparound=False -# cython: cdivision=True - -# Import the Python-level symbols of numpy -import numpy as np - -cimport numpy as np - -# Numpy must be initialized. When using numpy from C or Cython you must -# _always_ do that, or you will have segfaults -np.import_array() - -cimport cython -from libc.stdio cimport FILE, fclose, fopen -from libc.stdlib cimport free, malloc - - -cdef extern from "stlfile.h": - ctypedef unsigned int vertex_t - ctypedef unsigned int triangle_t - int loadstl(FILE *fp, char *comment, float **vertp, vertex_t *nvertp, vertex_t **trip, unsigned short **attrp, triangle_t *ntrip) - - -cdef class ArrayWrapper: - cdef void* data_ptr - cdef int size - cdef int dtype - - cdef set_data(self, int size, void* data_ptr, int dtype): - """ Constructor for the class. - Mallocs a memory buffer of size (n*sizeof(int)) and sets up - the numpy array. - Parameters: - ----------- - n -- Length of the array. - Data attributes: - ---------------- - data -- Pointer to an integer array. - alloc -- Size of the data buffer allocated. - """ - self.data_ptr = data_ptr - self.size = size - self.dtype = dtype - - def __array__(self): - cdef np.npy_intp shape[1] - shape[0] = self.size - ndarray = np.PyArray_SimpleNewFromData(1, shape, self.dtype, self.data_ptr) - return ndarray - - def __dealloc__(self): - """ Frees the array. """ - free(self.data_ptr) - - -def get_stl_data(str filename): - cdef: - FILE *fp - char comment[80] - float *vertp - vertex_t nverts - vertex_t *trip - unsigned short *attrp - triangle_t ntrip - np.ndarray verts_arr, tri_arr, attr_arr - int out - - fp = fopen(filename.encode('utf-8'), "rb") - if fp is NULL: - raise FileNotFoundError(2, "No such file or directory: '%s'" % filename) - - # Call the STL file loader - out = loadstl(fp, comment, &vertp, &nverts, &trip, &attrp, &ntrip) - fclose(fp) - if out == -2: - raise RuntimeError("Invalid or unrecognized STL file format.") - elif out != 0: - raise RuntimeError("Unable to load STL.") - - # wrap the vertex array - point_array_wrapper = ArrayWrapper() - point_array_wrapper.set_data(nverts*3, vertp, np.NPY_FLOAT32) - points = np.array(point_array_wrapper).reshape(nverts, 3) - - # wrap the vertex array - tri_array_wrapper = ArrayWrapper() - tri_array_wrapper.set_data(ntrip*3, trip, np.NPY_UINT32) - indices = np.array(tri_array_wrapper).reshape(ntrip, 3) - - return points, indices diff --git a/src/stl_reader/reader.py b/src/stl_reader/reader.py index 97b379b..1b62057 100644 --- a/src/stl_reader/reader.py +++ b/src/stl_reader/reader.py @@ -11,12 +11,6 @@ from pyvista.core.pointset import PolyData -def _check_stl_ascii(filename: str) -> bool: - """Check if a STL is ASCII.""" - with open(filename, "rb") as fid: - return fid.read(5) == b"solid" - - def _polydata_from_faces(points: npt.NDArray[float], faces: npt.NDArray[int]) -> "PolyData": """Generate a polydata from a faces array containing no padding and all triangles. @@ -99,9 +93,6 @@ def read(filename: str) -> Tuple[npt.NDArray[np.float32], npt.NDArray[np.uint32] [9005998, 9005999, 9005995]], dtype=uint32) """ - if _check_stl_ascii(filename): - raise RuntimeError("stl-reader only supports binary STL files") - return _stlfile_wrapper.get_stl_data(filename) diff --git a/src/stlfile.cpp b/src/stlfile.cpp index bc212a1..76fc278 100644 --- a/src/stlfile.cpp +++ b/src/stlfile.cpp @@ -23,7 +23,8 @@ THE SOFTWARE. Also modified by Alex Kaszynski 2023-2024: - Check is the file is ASCII and the elimination of stderr. -- Add nanobind interface +- Add nanobind interface. +- Incorporate an ASCII reader. */ @@ -83,6 +84,12 @@ static vertex_t vertex(uint32_t *verts, vertex_t nverts, vertex_t *vht, return ~(vertex_t)0; } +static uint32_t float_to_uint32(float value) { + uint32_t result; + std::memcpy(&result, &value, sizeof(float)); + return result; +} + STL_STATUS check_stl_format(FILE *fp) { if (fp == NULL) { printf("\n\tUnable to open the file"); @@ -104,10 +111,11 @@ STL_STATUS check_stl_format(FILE *fp) { if (strcmp(sixBytes, "solid ") == 0) { char line[100]; - fgets(line, 100, fp); - if (strncmp(line, "facet ", 6) == 0) { - return STL_ASCII; - } + // fgets(line, 100, fp); + // if (strncmp(line, "facet ", 6) == 0) { + // return STL_ASCII; + // } + return STL_ASCII; rewind(fp); } @@ -127,6 +135,328 @@ STL_STATUS check_stl_format(FILE *fp) { return STL_BINARY; } +// Fast string to float reader +static inline float fast_atof(const char *&p) { + bool neg = false; + if (*p == '-') { + neg = true; + ++p; + } else if (*p == '+') { + ++p; + } + + double integer_part = 0.0; + while (*p >= '0' && *p <= '9') { + integer_part = integer_part * 10.0 + (*p - '0'); + ++p; + } + + double fraction_part = 0.0; + double fraction_scale = 1.0; + if (*p == '.') { + ++p; + while (*p >= '0' && *p <= '9') { + fraction_part = fraction_part * 10.0 + (*p - '0'); + fraction_scale *= 10.0; + ++p; + } + } + + double exponent = 0.0; + if (*p == 'e' || *p == 'E') { + ++p; + bool exp_neg = false; + if (*p == '-') { + exp_neg = true; + ++p; + } else if (*p == '+') { + ++p; + } + while (*p >= '0' && *p <= '9') { + exponent = exponent * 10.0 + (*p - '0'); + ++p; + } + if (exp_neg) + exponent = -exponent; + } + + double result = (integer_part + fraction_part / fraction_scale); + if (exponent != 0.0) { + result *= pow(10.0, exponent); + } + + if (neg) + result = -result; + + return static_cast(result); +} + +int loadstl_ascii(FILE *fp, char *comment, float **vertp, vertex_t *nvertp, + vertex_t **trip, uint16_t **attrp, triangle_t *ntrip) { + // Read the entire file into memory + fseek(fp, 0, SEEK_END); + size_t fileSize = ftell(fp); + rewind(fp); + + char *fileBuffer = (char *)malloc(fileSize + 1); + if (!fileBuffer) { + fprintf(stderr, "Memory allocation failed for file buffer\n"); + return -1; + } + if (fread(fileBuffer, 1, fileSize, fp) != fileSize) { + fprintf(stderr, "Failed to read file into buffer\n"); + free(fileBuffer); + return -1; + } + fileBuffer[fileSize] = '\0'; // Null-terminate the buffer + + const char *ptr = fileBuffer; + const char *end = fileBuffer + fileSize; + + // Extract the comment from the first line + if (sscanf(ptr, "solid %79[^\n]", comment) != 1) { + strcpy(comment, ""); + } + // Move ptr to the end of the line + while (ptr < end && *ptr != '\n') + ptr++; + if (ptr < end) + ptr++; + + // First, estimate the number of triangles + triangle_t ntris_estimate = 0; + const char *scan_ptr = ptr; + while (scan_ptr < end) { + // Skip whitespace + while (scan_ptr < end && isspace(*scan_ptr)) + scan_ptr++; + + if (scan_ptr >= end) + break; + + if (strncmp(scan_ptr, "facet", 5) == 0) { + ntris_estimate++; + // Move to the end of the line + while (scan_ptr < end && *scan_ptr != '\n') + scan_ptr++; + if (scan_ptr < end) + scan_ptr++; + continue; + } + + // Move to the next line + while (scan_ptr < end && *scan_ptr != '\n') + scan_ptr++; + if (scan_ptr < end) + scan_ptr++; + } + + // Allocate initial memory for vertices and triangles + triangle_t ntris = 0; + vertex_t nverts = 0; + size_t tris_cap = ntris_estimate > 0 ? ntris_estimate : 1024; + vertex_t *tris = (vertex_t *)malloc(tris_cap * 3 * sizeof(vertex_t)); + size_t verts_cap = tris_cap * 3; + uint32_t *verts = (uint32_t *)malloc(verts_cap * 3 * sizeof(uint32_t)); + + // Set vhtcap based on estimated number of vertices + vertex_t vhtcap = nextpow2(verts_cap * 2); // Must be a power of two + vertex_t *vht = (vertex_t *)calloc(vhtcap, sizeof(vertex_t)); + if (!tris || !verts || !vht) { + fprintf(stderr, "Memory allocation failed\n"); + free(fileBuffer); + free(tris); + free(verts); + free(vht); + return -1; + } + + int v_idx = 0; + uint32_t v[3][3]; + uint32_t vert[3]; + + ptr = fileBuffer; + // Skip the first line (solid ...) + while (ptr < end && *ptr != '\n') + ptr++; + if (ptr < end) + ptr++; + + while (ptr < end) { + // Skip whitespace + while (ptr < end && isspace(*ptr)) + ptr++; + + if (ptr >= end) + break; + + if (strncmp(ptr, "facet", 5) == 0) { + ptr += 5; + // Skip to the end of the line + while (ptr < end && *ptr != '\n') + ptr++; + if (ptr < end) + ptr++; + v_idx = 0; + continue; + } + + if (strncmp(ptr, "vertex", 6) == 0) { + ptr += 6; + // Skip whitespace + while (ptr < end && isspace(*ptr)) + ptr++; + // Parse three floats + float x, y, z; + + x = fast_atof(ptr); + while (ptr < end && isspace(*ptr)) + ptr++; + + y = fast_atof(ptr); + while (ptr < end && isspace(*ptr)) + ptr++; + + z = fast_atof(ptr); + while (ptr < end && *ptr != '\n') + ptr++; + if (ptr < end) + ptr++; + + if (v_idx < 3) { + vert[0] = float_to_uint32(x); + vert[1] = float_to_uint32(y); + vert[2] = float_to_uint32(z); + copy96(v[v_idx], vert); + v_idx++; + } else { + fprintf(stderr, "Warning: more than 3 vertices in a facet\n"); + } + continue; + } + + if (strncmp(ptr, "endfacet", 8) == 0) { + ptr += 8; + // Skip to the end of the line + while (ptr < end && *ptr != '\n') + ptr++; + if (ptr < end) + ptr++; + if (v_idx == 3) { + vertex_t vi[3]; + for (int i = 0; i < 3; i++) { + vi[i] = vertex(verts, nverts, vht, vhtcap, v[i]); + if (vi[i] == ~(vertex_t)0) { + // Hash table is full, need to resize + vertex_t old_vhtcap = vhtcap; + vertex_t *old_vht = vht; + vhtcap *= 2; + vht = (vertex_t *)calloc(vhtcap, sizeof(vertex_t)); + if (!vht) { + fprintf(stderr, "Memory allocation failed for vht\n"); + free(fileBuffer); + free(tris); + free(verts); + free(old_vht); + return -1; + } + // Rehash existing entries + for (vertex_t idx = 0; idx < old_vhtcap; idx++) { + if (old_vht[idx] != 0) { + vertex_t vi_old = old_vht[idx] - 1; + uint32_t *vert_old = verts + 3 * vi_old; + vertex_t hash = final96(vert_old[0], vert_old[1], vert_old[2]); + vertex_t j; + for (j = 0; j < vhtcap; j++) { + vertex_t *vip_new = vht + ((hash + j) & (vhtcap - 1)); + if (*vip_new == 0) { + *vip_new = vi_old + 1; + break; + } + } + if (j == vhtcap) { + fprintf(stderr, "Failed to rehash during vht resize\n"); + free(fileBuffer); + free(tris); + free(verts); + free(old_vht); + free(vht); + return -1; + } + } + } + free(old_vht); + // Now retry inserting the current vertex + vi[i] = vertex(verts, nverts, vht, vhtcap, v[i]); + if (vi[i] == ~(vertex_t)0) { + fprintf(stderr, "Vertex hash table is full after resizing\n"); + free(fileBuffer); + free(tris); + free(verts); + free(vht); + return -1; + } + } + if (vi[i] == nverts) { + if (nverts >= verts_cap) { + verts_cap *= 2; + verts = + (uint32_t *)realloc(verts, verts_cap * 3 * sizeof(uint32_t)); + if (!verts) { + fprintf(stderr, "Failed to reallocate verts array\n"); + free(fileBuffer); + free(tris); + free(vht); + return -1; + } + } + copy96(verts + 3 * nverts, v[i]); + nverts++; + } + } + + if (ntris >= tris_cap) { + tris_cap *= 2; + tris = (vertex_t *)realloc(tris, tris_cap * 3 * sizeof(vertex_t)); + if (!tris) { + fprintf(stderr, "Failed to reallocate tris array\n"); + free(fileBuffer); + free(verts); + free(vht); + return -1; + } + } + tris[3 * ntris + 0] = vi[0]; + tris[3 * ntris + 1] = vi[1]; + tris[3 * ntris + 2] = vi[2]; + ntris++; + + v_idx = 0; + } else { + fprintf(stderr, "Incomplete facet, expected 3 vertices\n"); + } + continue; + } + + // Skip unrecognized lines + // Move to the end of the line + while (ptr < end && *ptr != '\n') + ptr++; + if (ptr < end) + ptr++; + } + + free(fileBuffer); + free(vht); + verts = (uint32_t *)realloc(verts, nverts * 3 * sizeof(uint32_t)); + *vertp = (float *)verts; + *nvertp = nverts; + *trip = tris; + *attrp = NULL; // ASCII STL does not have attributes + *ntrip = ntris; + return 0; +} + int loadstl(FILE *fp, char *comment, float **vertp, vertex_t *nvertp, vertex_t **trip, uint16_t **attrp, triangle_t *ntrip) { uint8_t buf[128]; @@ -143,6 +473,11 @@ int loadstl(FILE *fp, char *comment, float **vertp, vertex_t *nvertp, return -2; } + if (format_status == STL_ASCII) { + fseek(fp, 0, SEEK_SET); + return loadstl_ascii(fp, comment, vertp, nvertp, trip, attrp, ntrip); + } + // the comment and triangle count if (fread(buf, 84, 1, fp) != 1) { fprintf(stderr, "loadstl: short read at header\n"); diff --git a/tests/test_reader.py b/tests/test_reader.py index ffdc4b9..2765611 100644 --- a/tests/test_reader.py +++ b/tests/test_reader.py @@ -71,8 +71,7 @@ def test_read_binary() -> None: def test_read_ascii() -> None: - with pytest.raises(RuntimeError): - stl_reader.read(TEST_FILE_ASCII) + points, ind = stl_reader.read(TEST_FILE_ASCII) @pytest.mark.skipif(not PYVISTA_INSTALLED, reason="Requires PyVista") # type: ignore