MarchingNumPy

This package is a Python and NumPy implementation of the marching cubes and related algorithms.

Background

Data Visualisation

The visualisation of numeric data is a key component of computer graphics.

Data is often encountered as numeric values distibuted over a regular space. Sometimes these volumes data are called voxels. Visualisation of this numeric data is key to the interpretation and/or appreciation of the data.

Selected examples:

  • Medicine: medical imaging such as MRI and CT scanners.

  • Science: quantum mechanical calculations and data recording.

  • Cartography: linking data such as height, land usage, etc to map coordinates.

  • Entertainment: generation of artwork generated terrain and visualsation of metaballs.

An important method for the visualisation of volume data is to calculate an isosurface: a mesh surface that follows as closely as possible where the data crosses a certain level threshold.

Marching Cubes: Lorensen

This package includes an implementation of the Lorensen Marching Cubes algorithm using only python and numpy methods: MarchingCubesLorensen.marching_cubes_lorensen().

Lorensen and Cline described the Marching Cubes algorithm in 1987.

Marching Cubes joins all points where a volume of numeric values crosses a level threshold along each axis to create a triangular mesh. This algorithm considers each unit or cube of the data individually and performs the following operations:

  • Calculate Intersects.

    For each time that the volume value crosses the level threshold along each axis, interpolate the position of the intersection - this becomes a vertex for the output mesh.

  • Assign Types.

    Consider each corner of the volume unit and whether it is above or below the level threshold. Assign a unique identifier to each unit based upon the outcome of these logical binary tests. For a cube there are 8 vertices that gives \(2^8 = 1\newcommand*\ShiftLeft{\ll}8 = 256\) possible outcomes. However, these can be reduced to the same fundamental 14 types by symmetry.

  • Look Up Geometry.

    Use the calculated type to index a precalculated geometry table the defines how the calculated vertices should be joined using triangular geometry.

Drawbacks:

For intricate volume data, the simple reduction to 14 types of cube is insufficient to describe the isosurface. This ambiguity can lead to gaps within the mesh.

Alternative Marching Approaches:

  • Marching Tetrahedra

  • Marching Cubes 33

Marching Cubes in 2D: Marching Squares and Marching Triangles

This package includes an implementation of the Marching Squares and Marching Triangles algorithms using numpy: MarchingSquares.marching_squares() and MarchingTriangles.marching_triangles() / MarchingTriangles.marching_triangles_reversed().

The same logic can be applied 2-dimensionally to find isolines or contours for 2D values. The calculation of intersects is the same but now in 2D along two axes rather than 3D. The marching squares have 4 vertices that gives \(2^4 = 1\newcommand*\ShiftLeft{\ll}4 = 16\) possible types that are reduced to 7 by symmetry. It is trivial to draw up the look-up table by hand.

Marching Triangles is an exercise where the squares are split into triangles along their diagonals. Triangles can be split in either of two directions / or \. This approach results in a system with no ambiguity.

Marching NumPy

The importance of calculating isosurfaces this was means there are multiple implementations of the marching cubes algorithm. Additionally, the marching cubes algorithm is a common teaching tool in computer science.

  • The scikit-image package includes skimage.measure.marching_cubes() [API]

  • Within python, marching cubes is available from multiple sources: [PiPy “marching cubes”]

NumPy is a powerful package for manipulating data in python that is used widely for scientific computing and data analysis. NumPy provides a very efficient way to store arrayed data and a suite of methods to process that data. If NumPy methods are exploited correctly then calculations can be signficantly accelerated by exploting the underlying lower level functions. This package uses only NumPy methods to acheive “very fast” marching cubes.

References and Notes

NumPy
Publications
Wikipedia
Blogs and Tutorials
Other Marching Solutions

MarchingNumPy

MarchingNumPy.MarchingCubesLorensen.marching_cubes_lorensen(volume, level=0.0, *, interpolation='LINEAR', step_size=1, resolve_ambiguous=True)

Convert a 3d volume into an isosurface by splitting the volume into cubes and evaluating where the values in the volume cross a level threshold along the edges of the squares.

Ambiguity resolution is not available with this function.

Parameters:

  • volume (NDArray) – volume

  • level (float | NDArray) – level

Keyword Arguments:
Returns:

  • vertices (NDArray)

  • geometry (NDArray)

MarchingNumPy.MarchingSquares.marching_squares(volume, level=0.0, *, interpolation='LINEAR', step_size=1, resolve_ambiguous=True)

Convert a 2d volume into isolines by splitting the volume into squares and evaluating where the values in the volume cross a level threshold along the edges of the squares. Can perform a test where the volume type is ambiguous to resolve the ambiguity.

Parameters:

  • volume (NDArray) – volume

  • level (float | NDArray) – level

Keyword Arguments:

  • interpolation ({'LINEAR', 'HALFWAY', 'COSINE'}) – Interpolation method for FindIntersects.find_intersects().

  • step_size (int) – Defaults to 1.

  • ambiguity_resolution (bool) – Perform ambiguity resolution. Defaults to True.

Returns:

  • vertices (NDArray)

  • geometry (NDArray)

MarchingNumPy.MarchingTriangles.marching_triangles(volume, level=0.0, *, interpolation='LINEAR', step_size=1, resolve_ambiguous=True)

Convert a 2d volume into isolines by splitting the volume into triangles and evaluating where the values in the volume cross a level threshold along the edges of the triangles. The triangles are formed by splitting squares along the bottom left - top right diagonal. Ambiguity resolution is not relevant for this function.

Parameters:

  • volume (NDArray) – volume

  • level (float | NDArray) – level

Keyword Arguments:
Returns:

  • vertices (NDArray)

  • geometry (NDArray)

MarchingNumPy.MarchingTriangles.marching_triangles_reversed(volume, level=0.0, *, interpolation='LINEAR', step_size=1, resolve_ambiguous=True)

Convert a 2d volume into isolines by splitting the volume into triangles and evaluating where the values in the volume cross a level threshold along the edges of the triangles. The triangles are formed by splitting squares along the bottom right - top left diagonal. Ambiguity resolution is not relevant for this function.

Parameters:

  • volume (NDArray) – volume

  • level (float | NDArray) – level

Keyword Arguments:
Returns:

  • vertices (NDArray)

  • geometry (NDArray)

MarchingNumPy Algorithmic Steps

MarchingNumPy.Marching.marching_factory(*, nD, nEdges, intersect_slice_indexes, volume_type_slices, ambiguity_resolution=None, geometry_array, edge_direction, edge_delta)

Factory method for creating marching methods.

Keyword Arguments:
MarchingNumPy.FindIntersects.find_intersects(slice_indexes, volume, size_multiplier, volume_test=None, level=0.0, interpolation='LINEAR')

Calculate intersects in volume data.

This function considers values in a volume along slices as determined by slice_indexes. Each intersect where the volume crosses the level is recorded in the intersects NDArray. Each intersect has an associated unique id that is recorded in the intersect_ids NDArray.

The level is used to generate a zero-referenced volume and volume_test as below. If level is non-truthy (i.e. 0.0 or None) then volume is assumed to be zero-referenced. If volume_test is already calculated then this can be passed in and, as long as level is non-truthy, this part of the calculation will be skipped.

volume = volume - level
volume_test = volume >= 0

The crossing point is calculated according to value of interpolation:

“HALFWAY”

The crossing point is marked at 0.5 between the two corners. This is a fast interpolation but gives geometry with a block appearence.

“LINEAR”

The crossing point is interpolated using linear interpolation.

\[\begin{split}\begin{align*} y &= mx + c \\ y &= (v_{n+1}-v_{n})x + v_{n} \\ x &= \frac{y - v_{n}}{v_{n+1} - v_{n}} \end{align*}\end{split}\]

At the zeroed volume crossing point, \(y = 0\)

\[x = \frac{v_{n}}{v_{n} - v_{n+1}}\]
“COSINE”

The crossing point is interpolated using cosine interpolation.

\[\begin{split}\begin{align*} y &= \frac{h}{2}\cos{(\pi x)} + c \\ y &= \frac{v_{n}-v_{n+1}}{2} \cos({\pi x}) + \frac{v_{n} + v_{n+1}}{2} \\ x &= \frac{1}{\pi}\arccos\left({\frac{2y - v_{n} - v_{n+1}}{v_{n}-v_{n+1}}}\right) \end{align*}\end{split}\]

At the zeroed volume crossing point, \(y = 0\)

\[x = \frac{1}{\pi}\arccos\left({\frac{v_{n+1} + v_{n}}{v_{n+1}-v_{n}}}\right)\]
Parameters:

  • slice_indexes (Collection[Tuple[slice, ...], ...]) – A Collection of Tuple of Slices that describe where to look for the intersects.

  • volume (ArrayLike) – A n-dimensional array containing the volume data for testing. Passed through numpy.asarray().

  • size_multiplier (NDArray[Any]) – A multiplier to convert a coordinate and direction into a unique id.

  • volume_test (NDArray[numpy.bool_], default = None) – Optional bool-like results of the results of testing the volume against the level. The size must match the size of volume.

  • level (ArrayLike, default = 0.0) – The level against which volume is evaluated. Passed through numpy.asarray() if truthy.

  • interpolation ({"LINEAR", "HALFWAY", "COSINE"}) – Interpolation of the crossing point.

Returns:

  • intersects (NDArray[Intersect]) – A NDArray containing intersects, ordered by axis.

  • intersect_ids (NDArray[Intersect_id]) – A NDArray containing a unique id for each intersect in same order as intersects for indexing.

Raises:

ValueError – If the volume is zero-sized; or if the shape of the supplied volume_test is not consistent with volume; or if interpolation is an invalid value.

MarchingNumPy.VolumeTypes.volume_types(volume_test, slices, *, dtype=<class 'numpy.uint16'>)

Calculate the type of each unit of a volume based on the values at each corner.

The returned numpy array types is initalised to a numpy.zeros() array with shape n-1 for each n in volume_test shape. Then, types is updated depending on the values in volume_test (where volume_test is the result of an operation such as volume_test = volume >= level) in the direction stipulated by each of the slices. The function iterates through the collection of slices and types is updated with the index bit (shifted left every iteration of the loop).

Example slices for standard 3D and 2D volumes (with corresponding bit set as the comment). See IndexingTools.str_to_index() for a tool to generate these slices.

# Standard slices for a 3D volume
slices = [
    slice(None, -1), slice(None, -1), slice(None, -1), # 1:   x,   y,   z
    slice( 1, None), slice(None, -1), slice(None, -1), # 2:   x+1, y,   z
    slice( 1, None), slice( 1, None), slice(None, -1), # 4:   x+1, y+1, z
    slice(None, -1), slice( 1, None), slice(None, -1), # 8:   x,   y+1, z
    slice(None, -1), slice(None, -1), slice( 1, None), # 16:  x,   y,   z+1
    slice( 1, None), slice(None, -1), slice( 1, None), # 32:  x+1, y,   z+1
    slice( 1, None), slice( 1, None), slice( 1, None), # 64:  x+1, y+1, z+1
    slice(None, -1), slice( 1, None), slice( 1, None), # 128: x,   y+1, z+1
]

# Standard slices for a 2D volume
slices = [
    slice(None, -1), slice(None, -1), # 1: x,   y
    slice( 1, None), slice(None, -1), # 2: x+1, y
    slice( 1, None), slice( 1, None), # 4: x+1, y+1
    slice(None, -1), slice( 1, None), # 8: x,   y+1
]
Parameters:

  • volume_test (NDArray[numpy.bool_]) – A n-dimensional bool-like volume test. There must be at least one value in each dimension. Passed through numpy.asarray() with dtype=numpy.bool_.

  • slices (Collection[Tuple[slice, ...]]) – Collection of tuple of slices that is enumerated to give the type index and used to test where volume_test is True. See IndexingTools.str_to_index() for a tool to generate these slices.

Keyword Arguments:

dtype (Type[numpy.integer]) – The datatype used in the returned NDArray. Defaults to Types.VolumeType.

Returns:

types – A NDArray containing the results of the tests with of size with shape n-1 for n in each volume_test.shape.

Return type:

NDArray[numpy.integer]

Raises:

  • ValueError – If the arguments are not appropriate or the shapes of the supplied arrays are not consistent.

  • TypeError – If the supplied dtype is too small to hold the biggest possible type value.

MarchingNumPy.LookUpGeometry.look_up_geometry(types, geometry_array, edge_delta, edge_direction, size_multiplier, nDir=None)

Look up geometry from a look up table.

The value in types is used to lookup sets of edge number from geometry_array. The edge number information is used to lookup values in edge_direction and edge_delta. This can be reduced to an array of vertex_id_offsets once the size_multiplier is known. The current location is dot with the size_multiplier and added to vertex_id_offsets to calculate the vertex_id. The combined vertex_ids make up the geometry.

Parameters:

  • types (NDArray[numpy.unsignedinteger]) – A NDArray containing values to index the geometry_array. Passed through numpy.asarray().

  • geometry_array (NDArray[numpy.integer]) – The reference array containing the edge number information lookup values.

  • edge_delta (NDArray[numpy.integer]) – Edge delta coordinate lookup values.

  • edge_direction (NDArray[numpy.integer]) – Edge direction lookup values.

  • size_multiplier (NDArray[numpy.integer]) – Size multiplier to convert the edge number into an edge id.

  • nDir (int | None) – The number of vertex directions used to calculate the edge_id. Defaults to None when nDir = 1 + edge_id_direction.max() will be used.

Raises:

ValueError – If geometry_array has an incompatible shape.

Returns:

geometry – A NDArray of connected vertex_ids having shape: number of geometry found, number of vertices per geometry.

Return type:

NDArray[numpy.uint64]

MarchingNumPy.ConvertIndexes.convert_indexes(array, ordered_ids, *, method='NUMPY', **kwargs)

Replace values in array with the index of those values within an ordered_ids array.

Parameters:

  • array (ArrayLike) – array. Passed through numpy.asarray().

  • ordered_ids (ArrayLike) – ordered array. Passed through numpy.asarray().

  • method ({"NUMPY", "DICT"}) – use a numpy array for the lookup or use a dict for the lookup.

Returns:

convered array of values after looukp.

Return type:

NDArray

General Tools

MarchingNumPy.FindIntersects.vector_from_slices(from_slices, to_slices, absolute=False)

Calculates a direction vector based upon two Collections of slices.

Parameters:

  • from_slices (Collection[slice]) – from Collection of slices.

  • to_slices (Collection[slice]) – to Collection of slices.

  • absolute (bool) – If `True` returns absolute values. Default False.

Returns:

vector representing direction.

Return type:

NDArray[numpy.int8]

MarchingNumPy.ConvertIndexes.ndarray_dict_lookup(array, dictionary, *, dtype=None, **kwargs)

Vectorized lookup of array values within a dict.

Parameters:

  • array (NDArray) – array

  • dictionary (dict) – lookup values.

Keyword Arguments:

  • dtype (Type, optional) – dtype of the results. Defaults to None.

  • default (Any, optional) – the default value for dictionary.get().

Returns:

results of the dictionary lookup.

Return type:

NDArray

Raises:

KeyError – If a value in array is not in ordered_ids and default is not supplied.

MarchingNumPy.ConvertIndexes.ndarray_dict_ordered_lookup(array, ordered_ids, *, dtype=None, **kwargs)

Replace values in an array values with the index of those values within an ordered_ids array using ndarray_dict_lookup().

Parameters:

  • array (NDArray) – array.

  • ordered_ids (NDArray) – lookup values in order.

For keyword parameters, see ndarray_dict_lookup().

Returns:

results of the array lookup.

Return type:

NDArray

MarchingNumPy.ConvertIndexes.ndarray_numpy_ordered_lookup(array, ordered_ids, *, dtype=None)

Replace values in an array values with the index of those values within an ordered_ids array using :module:`numpy`.

N.B. Creates a large NDArray for the lookup!

Parameters:

  • array (NDArray) – array.

  • ordered_ids (NDArray) – lookup values in order.

Keyword Arguments:

dtype (Type, optional) – dtype of the results. Defaults to None.

Returns:

results of the array lookup.

Return type:

NDArray

Raises:

ValueError – If the array of values is too large for a numpy.uint64.

MarchingNumPy.Checking.assert_nd_array(a, nD, minsize=0)

Check the numpy array a has exactly nD dimensions and the shape in each dimension is greater than minsize.

Parameters:

  • a (ArrayLike) – A numpy ndarray. Passed through numpy.asarray().

  • nD (int) – The number of dimensions to test for.

  • minsize (int, default = 0) – The minimum required size in each dimension.

Raises:

ValueError – If any of the assertions fail.

MarchingNumPy.IndexingTools.int_or_none(i)

Convert to int, or on ValueError return None.

>>> print(int_or_none("3"))
3
>>> print(int_or_none(3.2))
3
>>> print(int_or_none(""))
None
>>> print(int_or_none("3.2"))
None
Parameters:

i (Any) – input for conversion.

Returns:

  • int – If the input can be conveted to an int.

  • None – If a ValueError arises during conversion.

MarchingNumPy.IndexingTools.str_to_index(string)

Convert a str representing an index to the approriate object in a similar way to the Python interpreter.

If the string contains a “,” then it will be interpreted as a multimensional NumPy-type index. It will be converted to a tuple of int, slices or Ellipsis using a single recursion of this function.

Otherwise the input wil be interpreted as a standard index. Single integer indexes are converted to int. “Ellipis” or “…” is converted to Ellipsis. If the string contains a “:” then it will be converted to a slice.

Square brackets are optional.

>>> str_to_index("1")
1
>>> str_to_index("[-1]")
-1
>>> str_to_index("[...]")
Ellipsis
>>> str_to_index("Ellipsis")
Ellipsis
>>> str_to_index("[1:]")
slice(1, None, None)
>>> str_to_index("[:-1]")
slice(None, -1, None)
>>> str_to_index("::2")
slice(None, None, 2)
>>> str_to_index("[1:, 1:, 1:]")
(slice(1, None, None), slice(1, None, None), slice(1, None, None))
>>> str_to_index("[..., 5]")
(Ellipsis, 5)
>>> str_to_index("[..., ::5]")
(Ellipsis, slice(None, None, 5))

Example to generate the slices for marching cubes VolumeTypes.volume_types().

list(
    map(
        str_to_index,
        [
            "[:-1, :-1, :-1]",  # x,   y  , z       1
            "[1: , :-1, :-1]",  # x+1, y  , z       2
            "[1: , 1: , :-1]",  # x+1, y+1, z       4
            "[:-1, 1: , :-1]",  # x,   y+1, z       8
            "[:-1, :-1, 1: ]",  # x,   y  , z+1    16
            "[1: , :-1, 1: ]",  # x+1, y  , z+1    32
            "[1: , 1: , 1: ]",  # x+1, y+1, z+1    64
            "[:-1, 1: , 1: ]",  # x,   y+1, z+1   128
        ],
    )
)
Parameters:

string (str) – input string (square brackets are optional).

Returns:

  • tuple[int, Ellipsis, slice] – if “,” in the input then it is interpreted as a multidimensional index.

  • int – if the input can be converted to an int.

  • Ellipsis – if the input is “…” or “Ellipsis”.

  • slice – if “:” in the input then it will be converted to a slice.

MarchingNumPy.IndexingTools.str_to_index_tuple(*string)

See: str_to_index() Takes one or more str values and returns a tuple of the result of mapping str_to_index() to each.

MarchingCuPy

MarchingCuPy is a direct conversion of MarchingNumPy to use CuPy instead of NumPy. The speedup by using CuPy where available is incredible - although there is a cost associated with inital setup.

The module code is created automatically by adding import cupy and replacing numpy. array creation routines with their cupy. equivalents.

Substitutions:

  • numpy.asarray → cupy.asarray

  • numpy.concatenate → cupy.concatenate

  • numpy.full_like → cupy.full_like

  • numpy.nonzero → cupy.nonzero

  • numpy.ones → cupy.ones

  • numpy.zeros → cupy.zeros

The majority of operations are methods on array objects so are converted implicitly. In other cases the cupy. attribute is an alias for the numpy. attribute.

The code in MarchingNumPy is idential to that in MarchingCuPy with a few subtle differences in operation.

Operational differences:

  • cupy does not support the where= keyword argument in methods so an alternative is used

  • cupy.einsum() is exceptionally slow so an alternative is used (also cupy.dot())

Indices and tables