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:

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 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:

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 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:

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 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.