ilastik/marching_cubes: A Small, Fast C++ Iso-Surface Extractor for Python

What It Does and Who It Is For

This project is a standalone implementation of the marching cubes algorithm, written in C++ with Python bindings. It creates a 3D iso-surface from a 3D volume, which is a mesh of triangles that represents a surface where the volume data crosses a given threshold. The library was originally developed for the volumina 3D viewer, so its primary audience is developers who need to visualize volumetric data, such as medical images or scientific simulations, within a Python application. It is not a general-purpose mesh processing library; it does one thing: extract a surface mesh from a scalar volume. If you need to render a surface from a CT scan or a fluid simulation, this library can do that with minimal dependencies. If you need more complex operations like decimation or boolean operations, you will need to look elsewhere.

How It Works: The Marching Cubes Algorithm and Smoothing

The core mechanism is the classic marching cubes algorithm. The library takes a 3D array of values and a threshold. It then processes each cube of eight neighboring voxels, determines which edges the iso-surface crosses, and generates triangles accordingly. The README shows that the `march` function returns three arrays: vertices, normals, and faces. The normals are computed as part of the extraction, which is a key point because some implementations require you to compute normals separately. The library also includes a smoothing step, controlled by an integer argument. The README notes that the number of smoothing rounds was off by one in version 0.2, and that in version 0.3 they fixed it. So `march(volume, 3)` performs three smoothing rounds, but in older versions you would have needed `march(volume, 4)` to get the same result. This smoothing is likely a Laplacian-like smoothing that moves vertices to reduce mesh noise, but the exact algorithm is not documented in the provided material. The output is a set of vertices and faces that can be directly fed into a rendering library like pyqtgraph.

Getting It Running: Installation and Basic Usage

The README does not provide explicit installation instructions, but since it is a Python package with C++ bindings, you would typically install it via `pip install marching_cubes` or build from source. The example usage is clear: you import `march` from `marching_cubes`, load a 3D numpy array, and call `march(volume, isovalue)`. The isovalue is the threshold; the README says 'extract the mesh where the values are larger than or equal to 1', so the second argument is the threshold, not the number of smoothing rounds. Wait, the example shows `march(volume, 0)` and comments 'zero smoothing rounds'. That is confusing. Looking closer, the signature is `march(volume, smoothing_rounds)`, and the threshold is presumably set inside the volume? Actually, the README says 'extract the mesh where the values are larger than or equal to 1', but the call is `march(volume, 0)`. That suggests the second argument is the smoothing rounds, and the threshold is hardcoded or determined by the data? No, that cannot be right. The text says 'extract the mesh where the values are larger than or equal to 1', but then 'march(volume, 0)'. It is possible that the threshold is implicit in the volume, i.e., the volume is already binary or the isovalue is 0? Actually, the volume is uint8, and the threshold might be 0.5? I cannot confirm. The README is ambiguous. I will state that the example uses `march(volume, 0)` for zero smoothing rounds, and that the threshold is presumably derived from the data or a separate parameter, but the documentation does not clarify. This is a genuine limitation in the documentation.

Memory Layout and the Fortran Order Fix

One of the most important changes in version 0.3 was fixing the expectation of Fortran order data. The release notes say that previously, the library expected data in Fortran order, and a workaround was to supply `volume.T` with C-order data. As of version 0.3, `march` works with both C-order and F-order data. This is a significant usability improvement because NumPy arrays are C-order by default, and requiring Fortran order would be a common pitfall. However, the note also warns that this change breaks compatibility with older versions: if you were using `volume.T` as a workaround, you should stop doing that. This is a concrete example of how the library's behavior changed and how you must update your code when upgrading. For a library that is meant to be lightweight, this kind of fix reduces friction, but it also means that you need to check which version you are using and adjust your code accordingly.

Performance Claims and What They Mean

The README includes a performance example on a 128x128x128 volume. It reports that extracting a mesh with 82,464 vertices and 165,048 triangles took 0.158 seconds, with an additional 0.254 seconds for three smoothing rounds, on an AMD A8-4500M 1.9 GHz CPU. That processor is from 2012, so on modern hardware the times would be lower. These numbers are specific to that volume size and hardware, so they are not a benchmark you can extrapolate to other scenarios. But they do give a rough idea: the algorithm is fast enough for interactive use on moderate volumes. The README also compares against scikit-image, claiming it is 10 to 20 times slower and produces incorrect normals. That comparison is from 2016, so it may be outdated. Still, the claim suggests that for performance-critical applications, a C++ implementation like this one has an advantage over pure Python or Cython implementations. You should not take these numbers as gospel, but they indicate the library was designed with speed in mind.

Why Not skimage or VTK: The Original Trade-offs

The README includes a section titled 'Why not these?' that explains why the authors chose to write their own implementation instead of using existing libraries. For scikit-image, they cite two reasons: it is too slow (factor 10-20) and it produces incorrect normals, leading to ugly shading. For VTK, they cite a huge dependency and lack of Python 3 support at the time of writing. They note that VTK has since added Python 3 support, so that objection is less relevant today. This section is useful because it frames the library's design goals: speed and correct normals, with minimal dependencies. If you are considering using this library, you should evaluate whether these trade-offs still hold for your use case. For example, if you are already using VTK for other purposes, adding it for mesh extraction might be acceptable, and you might not need this library. On the other hand, if you want a small, self-contained solution, this library fits.

Limitations and Known Issues

The README explicitly states that this library is 'work in progress'. The only listed todo is to 'repair workaround for volumes with different shapes'. That suggests there may be issues with volumes that are not cubical or have non-uniform dimensions. The example uses a 128x128x128 volume, so it is untested for rectangular volumes. Also, the threshold parameter is not clearly documented. The example says 'extract the mesh where the values are larger than or equal to 1', but the call is `march(volume, 0)`. That mismatch is confusing. I cannot tell from the material whether the second argument is the isovalue or the smoothing rounds. The README says 'zero smoothing rounds' for `march(volume, 0)`, so it is likely that the second argument is smoothing rounds, and the threshold is either hardcoded or determined by the volume's data type. That would be a limitation: you cannot specify an arbitrary threshold. I must state this ambiguity explicitly. Another limitation is that the library only outputs vertices, normals, and faces; it does not provide any mesh simplification or quality control. If you need to reduce triangle count, you would have to do that separately.

Maintenance and License Considerations

The project is licensed under BSD-3-Clause, which is permissive and allows commercial use with attribution. The last push was on 2025-07-25, which is recent, and there is a maintenance release 0.3.2 from the same date. The release history shows a gap of about five years between 0.3 and 0.3.1, but then a quick follow-up in 2025. That pattern suggests that the library is maintained on an as-needed basis, likely tied to the needs of the ilastik project. The README does not include a contribution guide or a list of maintainers, so you should not expect frequent updates. The maintenance cost for you is low: the API is small and stable, and the main risk is that the library does not support newer Python versions or numpy changes. You should verify that the Python bindings work with your Python version, as the README does not mention Python version compatibility. The BSD-3-Clause license is favorable, but you should read the full license text, which is not provided here.

Editorial conclusion

Adopt ilastik/marching_cubes if you need a lightweight, fast iso-surface extractor for 3D volumes in Python and you can accept its narrow scope: no advanced features like surface normals from gradients or multi-material extraction. Do not use it if you require a full-featured visualization pipeline or if your volumes are large enough that a few hundred milliseconds per mesh matters and you can afford a heavier dependency like VTK. Before adopting, verify that your data is C-order or F-order and that you understand the smoothing round semantics: version 0.3 and later expect C-order data, and the number of smoothing rounds is off by one compared to pre-0.2 behavior. Also confirm that the mesh output is suitable for your rendering target, since the normals are computed internally and may not match what you expect from gradient-based methods.

Frequently asked questions

Official sources

1. Official README
2. Project repository
3. Release notes