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 alevel
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:
interpolation ({'LINEAR', 'HALFWAY', 'COSINE'}) – Interpolation method for
FindIntersects.find_intersects()
.step_size (int) – Defaults to 1.
- 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 alevel
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 alevel
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:
interpolation ({'LINEAR', 'HALFWAY', 'COSINE'}) – Interpolation method for
FindIntersects.find_intersects()
.step_size (int) – Defaults to 1.
- 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 alevel
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:
interpolation ({'LINEAR', 'HALFWAY', 'COSINE'}) – Interpolation method for
FindIntersects.find_intersects()
.step_size (int) – Step size. Defaults to 1.
- 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:
nD (int) – Number of dimensions.
nDir (int) – Number of edge directions.
intersection_slice_indexes (Collection[Any]) – Intersection slice indexes for
FindIntersects.find_intersects()
.volume_type_slices (Collection[Tuple[slice, ...]]) – Slices for
VolumeTypes.volume_types()
.geometry_array (NDArray[numpy.integer[Any]]) – Geometry lookup array for
LookUpGeometry.look_up_geometry()
.edge_id_direction (NDArray[numpy.unsignedinteger[Any]]) – Edge ID direction for
LookUpGeometry.look_up_geometry()
.edge_id_delta (NDArray[numpy.unsignedinteger[Any]]) – Edge ID delta for
LookUpGeometry.look_up_geometry()
.ambiguity_reaolution (Optional[Callable[[NDArray, NDArray], NDArray["bool"]]]) – Callable function to modify types in ambiguous sitiations. Defaults to None.
- 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 ofvolume
.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 suppliedvolume_test
is not consistent withvolume
; or ifinterpolation
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 anumpy.zeros()
array with shapen-1
for eachn
involume_test
shape. Then,types
is updated depending on the values involume_test
(wherevolume_test
is the result of an operation such asvolume_test = volume >= level
) in the direction stipulated by each of theslices
. The function iterates through the collection ofslices
andtypes
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). SeeIndexingTools.str_to_index()
for a tool to generate theseslices
.# 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
isTrue
. SeeIndexingTools.str_to_index()
for a tool to generate theseslices
.
- 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
forn
in eachvolume_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 fromgeometry_array
. The edge number information is used to lookup values inedge_direction
andedge_delta
. This can be reduced to an array of vertex_id_offsets once the size_multiplier is known. The current location is dot with thesize_multiplier
and added tovertex_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 throughnumpy.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 whennDir = 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 anordered_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 inordered_ids
anddefault
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 returnNone
.>>> 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
ofint
,slices
orEllipsis
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 toEllipsis
. If the string contains a “:” then it will be converted to aslice
.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 mappingstr_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 replacingnumpy.
array creation routines with theircupy.
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 thenumpy.
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 (alsocupy.dot()
)