mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-15 11:25:53 +08:00
Merge pull request #469 from JDWarner/add_marching_cubes
3D Marching Cubes
This commit is contained in:
Emmanuelle Gouillart
+1
-1
@@ -114,7 +114,7 @@
|
||||
|
||||
- Joshua Warner
|
||||
Multichannel random walker segmentation, unified peak finder backend,
|
||||
n-dimensional array padding, bug and doc fixes.
|
||||
n-dimensional array padding, marching cubes, bug and doc fixes.
|
||||
|
||||
- Petter Strandmark
|
||||
Perimeter calculation in regionprops.
|
||||
|
||||
@@ -51,6 +51,9 @@ Library:
|
||||
Extension: skimage.measure._moments
|
||||
Sources:
|
||||
skimage/measure/_moments.pyx
|
||||
Extension: skimage.measure._marching_cubes_cy
|
||||
Sources:
|
||||
skimage/measure/_marching_cubes_cy.pyx
|
||||
Extension: skimage.graph._mcp
|
||||
Sources:
|
||||
skimage/graph/_mcp.pyx
|
||||
|
||||
@@ -0,0 +1,56 @@
|
||||
"""
|
||||
==============
|
||||
Marching Cubes
|
||||
==============
|
||||
|
||||
Marching cubes is an algorithm to extract a 2D surface mesh from a 3D volume.
|
||||
This can be conceptualized as a 3D generalization of isolines on topographical
|
||||
or weather maps. It works by iterating across the volume, looking for regions
|
||||
which cross the level of interest. If such regions are found, triangulations
|
||||
are generated and added to an output mesh. The final result is a set of
|
||||
vertices and a set of triangular faces.
|
||||
|
||||
The algorithm requires a data volume and an isosurface value. For example, in
|
||||
CT imaging Hounsfield units of +700 to +3000 represent bone. So, one potential
|
||||
input would be a reconstructed CT set of data and the value +700, to extract
|
||||
a mesh for regions of bone or bone-like density.
|
||||
|
||||
This implementation also works correctly on anisotropic datasets, where the
|
||||
voxel spacing is not equal for every spatial dimension, through use of the
|
||||
`sampling` kwarg.
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.mplot3d import Axes3D
|
||||
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
|
||||
|
||||
from skimage import measure
|
||||
from skimage.draw import ellipsoid
|
||||
|
||||
# Generate a level set about zero of two identical ellipsoids in 3D
|
||||
ellip_base = ellipsoid(6, 10, 16, levelset=True)
|
||||
ellip_double = np.concatenate((ellip_base[:-1, ...],
|
||||
ellip_base[2:, ...]), axis=0)
|
||||
|
||||
# Use marching cubes to obtain the surface mesh of these ellipsoids
|
||||
verts, faces = measure.marching_cubes(ellip_double, 0)
|
||||
|
||||
# Display resulting triangular mesh using Matplotlib. This can also be done
|
||||
# with mayavi (see skimage.measure.marching_cubes docstring).
|
||||
fig = plt.figure(figsize=(10, 12))
|
||||
ax = fig.add_subplot(111, projection='3d')
|
||||
|
||||
# Fancy indexing: `verts[faces]` to generate a collection of triangles
|
||||
mesh = Poly3DCollection(verts[faces])
|
||||
ax.add_collection3d(mesh)
|
||||
|
||||
ax.set_xlabel("x-axis: a = 6 per ellipsoid")
|
||||
ax.set_ylabel("y-axis: b = 10")
|
||||
ax.set_zlabel("z-axis: c = 16")
|
||||
|
||||
ax.set_xlim(0, 24) # a = 6 (times two for 2nd ellipsoid)
|
||||
ax.set_ylim(0, 20) # b = 10
|
||||
ax.set_zlim(0, 32) # c = 16
|
||||
|
||||
plt.show()
|
||||
@@ -1,11 +1,14 @@
|
||||
from .draw import circle, ellipse, set_color
|
||||
from ._draw import line, polygon, ellipse_perimeter, circle_perimeter, \
|
||||
bezier_segment
|
||||
from .draw3d import ellipsoid, ellipsoid_stats
|
||||
|
||||
__all__ = ['line',
|
||||
'polygon',
|
||||
'ellipse',
|
||||
'ellipse_perimeter',
|
||||
'ellipsoid',
|
||||
'ellipsoid_stats',
|
||||
'circle',
|
||||
'circle_perimeter',
|
||||
'set_color']
|
||||
|
||||
@@ -0,0 +1,117 @@
|
||||
# coding: utf-8
|
||||
import numpy as np
|
||||
from scipy.special import (ellipkinc as ellip_F, ellipeinc as ellip_E)
|
||||
|
||||
|
||||
def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False):
|
||||
"""
|
||||
Generates ellipsoid with semimajor axes aligned with grid dimensions
|
||||
on grid with specified `sampling`.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
a : float
|
||||
Length of semimajor axis aligned with x-axis.
|
||||
b : float
|
||||
Length of semimajor axis aligned with y-axis.
|
||||
c : float
|
||||
Length of semimajor axis aligned with z-axis.
|
||||
sampling : tuple of floats, length 3
|
||||
Sampling in (x, y, z) spatial dimensions.
|
||||
levelset : bool
|
||||
If True, returns the level set for this ellipsoid (signed level
|
||||
set about zero, with positive denoting interior) as np.float64.
|
||||
False returns a binarized version of said level set.
|
||||
|
||||
Returns
|
||||
-------
|
||||
ellip : (N, M, P) array
|
||||
Ellipsoid centered in a correctly sized array for given `sampling`.
|
||||
Boolean dtype unless `levelset=True`, in which case a float array is
|
||||
returned with the level set above 0.0 representing the ellipsoid.
|
||||
|
||||
"""
|
||||
if (a <= 0) or (b <= 0) or (c <= 0):
|
||||
raise ValueError('Parameters a, b, and c must all be > 0')
|
||||
|
||||
offset = np.r_[1, 1, 1] * np.r_[sampling]
|
||||
|
||||
# Calculate limits, and ensure output volume is odd & symmetric
|
||||
low = np.ceil((- np.r_[a, b, c] - offset))
|
||||
high = np.floor((np.r_[a, b, c] + offset + 1))
|
||||
|
||||
for dim in range(3):
|
||||
if (high[dim] - low[dim]) % 2 == 0:
|
||||
low[dim] -= 1
|
||||
num = np.arange(low[dim], high[dim], sampling[dim])
|
||||
if 0 not in num:
|
||||
low[dim] -= np.max(num[num < 0])
|
||||
|
||||
# Generate (anisotropic) spatial grid
|
||||
x, y, z = np.mgrid[low[0]:high[0]:sampling[0],
|
||||
low[1]:high[1]:sampling[1],
|
||||
low[2]:high[2]:sampling[2]]
|
||||
|
||||
if not levelset:
|
||||
arr = ((x / float(a)) ** 2 +
|
||||
(y / float(b)) ** 2 +
|
||||
(z / float(c)) ** 2) <= 1
|
||||
else:
|
||||
arr = ((x / float(a)) ** 2 +
|
||||
(y / float(b)) ** 2 +
|
||||
(z / float(c)) ** 2) - 1
|
||||
|
||||
return arr
|
||||
|
||||
|
||||
def ellipsoid_stats(a, b, c, sampling=(1., 1., 1.)):
|
||||
"""
|
||||
Calculates analytical surface area and volume for ellipsoid with
|
||||
semimajor axes aligned with grid dimensions of specified `sampling`.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
a : float
|
||||
Length of semimajor axis aligned with x-axis.
|
||||
b : float
|
||||
Length of semimajor axis aligned with y-axis.
|
||||
c : float
|
||||
Length of semimajor axis aligned with z-axis.
|
||||
sampling : tuple of floats, length 3
|
||||
Sampling in (x, y, z) spatial dimensions.
|
||||
|
||||
Returns
|
||||
-------
|
||||
vol : float
|
||||
Calculated volume of ellipsoid.
|
||||
surf : float
|
||||
Calculated surface area of ellipsoid.
|
||||
|
||||
"""
|
||||
if (a <= 0) or (b <= 0) or (c <= 0):
|
||||
raise ValueError('Parameters a, b, and c must all be > 0')
|
||||
|
||||
# Calculate volume & surface area
|
||||
# Surface calculation requires a >= b >= c and a != c.
|
||||
abc = [a, b, c]
|
||||
abc.sort(reverse=True)
|
||||
a = abc[0]
|
||||
b = abc[1]
|
||||
c = abc[2]
|
||||
|
||||
# Volume
|
||||
vol = 4 / 3. * np.pi * a * b * c
|
||||
|
||||
# Analytical ellipsoid surface area
|
||||
phi = np.arcsin((1. - (c ** 2 / (a ** 2.))) ** 0.5)
|
||||
d = float((a ** 2 - c ** 2) ** 0.5)
|
||||
m = (a ** 2 * (b ** 2 - c ** 2) /
|
||||
float(b ** 2 * (a ** 2 - c ** 2)))
|
||||
F = ellip_F(phi, m)
|
||||
E = ellip_E(phi, m)
|
||||
|
||||
surf = 2 * np.pi * (c ** 2 +
|
||||
b * c ** 2 / d * F +
|
||||
b * d * E)
|
||||
|
||||
return vol, surf
|
||||
@@ -0,0 +1,104 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_equal, assert_allclose
|
||||
|
||||
from skimage.draw import ellipsoid, ellipsoid_stats
|
||||
|
||||
|
||||
def test_ellipsoid_bool():
|
||||
test = ellipsoid(2, 2, 2)[1:-1, 1:-1, 1:-1]
|
||||
test_anisotropic = ellipsoid(2, 2, 4, sampling=(1., 1., 2.))
|
||||
test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]
|
||||
|
||||
expected = np.array([[[0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0],
|
||||
[0, 0, 1, 0, 0],
|
||||
[0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0]],
|
||||
|
||||
[[0, 0, 0, 0, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 0]],
|
||||
|
||||
[[0, 0, 1, 0, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[1, 1, 1, 1, 1],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 0, 1, 0, 0]],
|
||||
|
||||
[[0, 0, 0, 0, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 1, 1, 1, 0],
|
||||
[0, 0, 0, 0, 0]],
|
||||
|
||||
[[0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0],
|
||||
[0, 0, 1, 0, 0],
|
||||
[0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0]]])
|
||||
|
||||
assert_array_equal(test, expected.astype(bool))
|
||||
assert_array_equal(test_anisotropic, expected.astype(bool))
|
||||
|
||||
|
||||
def test_ellipsoid_levelset():
|
||||
test = ellipsoid(2, 2, 2, levelset=True)[1:-1, 1:-1, 1:-1]
|
||||
test_anisotropic = ellipsoid(2, 2, 4, sampling=(1., 1., 2.),
|
||||
levelset=True)
|
||||
test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]
|
||||
|
||||
expected = np.array([[[ 2. , 1.25, 1. , 1.25, 2. ],
|
||||
[ 1.25, 0.5 , 0.25, 0.5 , 1.25],
|
||||
[ 1. , 0.25, 0. , 0.25, 1. ],
|
||||
[ 1.25, 0.5 , 0.25, 0.5 , 1.25],
|
||||
[ 2. , 1.25, 1. , 1.25, 2. ]],
|
||||
|
||||
[[ 1.25, 0.5 , 0.25, 0.5 , 1.25],
|
||||
[ 0.5 , -0.25, -0.5 , -0.25, 0.5 ],
|
||||
[ 0.25, -0.5 , -0.75, -0.5 , 0.25],
|
||||
[ 0.5 , -0.25, -0.5 , -0.25, 0.5 ],
|
||||
[ 1.25, 0.5 , 0.25, 0.5 , 1.25]],
|
||||
|
||||
[[ 1. , 0.25, 0. , 0.25, 1. ],
|
||||
[ 0.25, -0.5 , -0.75, -0.5 , 0.25],
|
||||
[ 0. , -0.75, -1. , -0.75, 0. ],
|
||||
[ 0.25, -0.5 , -0.75, -0.5 , 0.25],
|
||||
[ 1. , 0.25, 0. , 0.25, 1. ]],
|
||||
|
||||
[[ 1.25, 0.5 , 0.25, 0.5 , 1.25],
|
||||
[ 0.5 , -0.25, -0.5 , -0.25, 0.5 ],
|
||||
[ 0.25, -0.5 , -0.75, -0.5 , 0.25],
|
||||
[ 0.5 , -0.25, -0.5 , -0.25, 0.5 ],
|
||||
[ 1.25, 0.5 , 0.25, 0.5 , 1.25]],
|
||||
|
||||
[[ 2. , 1.25, 1. , 1.25, 2. ],
|
||||
[ 1.25, 0.5 , 0.25, 0.5 , 1.25],
|
||||
[ 1. , 0.25, 0. , 0.25, 1. ],
|
||||
[ 1.25, 0.5 , 0.25, 0.5 , 1.25],
|
||||
[ 2. , 1.25, 1. , 1.25, 2. ]]])
|
||||
|
||||
assert_allclose(test, expected)
|
||||
assert_allclose(test_anisotropic, expected)
|
||||
|
||||
|
||||
def test_ellipsoid_stats():
|
||||
# Test comparison values generated by Wolfram Alpha
|
||||
vol, surf = ellipsoid_stats(6, 10, 16)
|
||||
assert(round(1280 * np.pi, 4) == round(vol, 4))
|
||||
assert(1383.28 == round(surf, 2))
|
||||
|
||||
# Test when a <= b <= c does not hold
|
||||
vol, surf = ellipsoid_stats(16, 6, 10)
|
||||
assert(round(1280 * np.pi, 4) == round(vol, 4))
|
||||
assert(1383.28 == round(surf, 2))
|
||||
|
||||
# Larger test to ensure reliability over broad range
|
||||
vol, surf = ellipsoid_stats(17, 27, 169)
|
||||
assert(round(103428 * np.pi, 4) == round(vol, 4))
|
||||
assert(37426.3 == round(surf, 1))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
np.testing.run_module_suite()
|
||||
@@ -1,4 +1,5 @@
|
||||
from .find_contours import find_contours
|
||||
from ._marching_cubes import marching_cubes, mesh_surface_area
|
||||
from ._regionprops import regionprops, perimeter
|
||||
from ._structural_similarity import structural_similarity
|
||||
from ._polygon import approximate_polygon, subdivide_polygon
|
||||
@@ -21,4 +22,7 @@ __all__ = ['find_contours',
|
||||
'moments',
|
||||
'moments_central',
|
||||
'moments_normalized',
|
||||
'moments_hu']
|
||||
'moments_hu',
|
||||
'sum_blocks',
|
||||
'marching_cubes',
|
||||
'mesh_surface_area']
|
||||
|
||||
@@ -0,0 +1,157 @@
|
||||
import numpy as np
|
||||
from . import _marching_cubes_cy
|
||||
|
||||
|
||||
def marching_cubes(volume, level, sampling=(1., 1., 1.)):
|
||||
"""
|
||||
Marching cubes algorithm to find iso-valued surfaces in 3d volumetric data
|
||||
|
||||
Parameters
|
||||
----------
|
||||
volume : (M, N, P) array of doubles
|
||||
Input data volume to find isosurfaces. Will be cast to `np.float64`.
|
||||
level : float
|
||||
Contour value to search for isosurfaces in `volume`.
|
||||
sampling : length-3 tuple of floats
|
||||
Voxel spacing in spatial dimensions corresponding to numpy array
|
||||
indexing dimensions (M, N, P) as in `volume`.
|
||||
|
||||
Returns
|
||||
-------
|
||||
verts : (V, 3) array
|
||||
Spatial coordinates for V unique mesh vertices. Coordinate order
|
||||
matches input `volume` (M, N, P).
|
||||
faces : (F, 3) array
|
||||
Define triangular faces via referencing vertex indices from ``verts``.
|
||||
This algorithm specifically outputs triangles, so each face has
|
||||
exactly three indices.
|
||||
|
||||
Notes
|
||||
-----
|
||||
The marching cubes algorithm is implemented as described in [1]_.
|
||||
A simple explanation is available here::
|
||||
|
||||
http://www.essi.fr/~lingrand/MarchingCubes/algo.html
|
||||
|
||||
There are several known ambiguous cases in the marching cubes algorithm.
|
||||
Using point labeling as in [1]_, Figure 4, as shown:
|
||||
|
||||
v8 ------ v7
|
||||
/ | / | y
|
||||
/ | / | ^ z
|
||||
v4 ------ v3 | | /
|
||||
| v5 ----|- v6 |/ (note: NOT right handed!)
|
||||
| / | / ----> x
|
||||
| / | /
|
||||
v1 ------ v2
|
||||
|
||||
Most notably, if v4, v8, v2, and v6 are all >= `level` (or any
|
||||
generalization of this case) two parallel planes are generated by this
|
||||
algorithm, separating v4 and v8 from v2 and v6. An equally valid
|
||||
interpretation would be a single connected thin surface enclosing all
|
||||
four points. This is the best known ambiguity, though there are others.
|
||||
|
||||
This algorithm does not attempt to resolve such ambiguities; it is a naive
|
||||
implementation of marching cubes as in [1]_, but may be a good beginning
|
||||
for work with more recent techniques (Dual Marching Cubes, Extended
|
||||
Marching Cubes, Cubic Marching Squares, etc.).
|
||||
|
||||
Because of interactions between neighboring cubes, the isosurface(s)
|
||||
generated by this algorithm are NOT guaranteed to be closed, particularly
|
||||
for complicated contours. Furthermore, this algorithm does not guarantee
|
||||
a single contour will be returned. Indeed, ALL isosurfaces which cross
|
||||
`level` will be found, regardless of connectivity.
|
||||
|
||||
The output is a triangular mesh consisting of a set of unique vertices and
|
||||
connecting triangles. The order of these vertices and triangles in the
|
||||
output list is determined by the position of the smallest ``x,y,z`` (in
|
||||
lexicographical order) coordinate in the contour. This is a side-effect
|
||||
of how the input array is traversed, but can be relied upon.
|
||||
|
||||
To quantify the area of an isosurface generated by this algorithm, pass
|
||||
the outputs directly into `skimage.measure.mesh_surface_area`.
|
||||
|
||||
Regarding visualization of algorithm output, the ``mayavi`` package
|
||||
is recommended. To contour a volume named `myvolume` about the level 0.0:
|
||||
|
||||
>>> from mayavi import mlab
|
||||
>>> verts, tris = marching_cubes(myvolume, 0.0, (1., 1., 2.))
|
||||
>>> mlab.triangular_mesh([vert[0] for vert in verts],
|
||||
[vert[1] for vert in verts],
|
||||
[vert[2] for vert in verts],
|
||||
tris)
|
||||
>>> mlab.show()
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Lorensen, William and Harvey E. Cline. Marching Cubes: A High
|
||||
Resolution 3D Surface Construction Algorithm. Computer Graphics
|
||||
(SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170).
|
||||
|
||||
See Also
|
||||
--------
|
||||
skimage.measure.mesh_surface_area
|
||||
|
||||
"""
|
||||
# Check inputs and ensure `volume` is C-contiguous for memoryviews
|
||||
if volume.ndim != 3:
|
||||
raise ValueError("Input volume must have 3 dimensions.")
|
||||
if level < volume.min() or level > volume.max():
|
||||
raise ValueError("Contour level must be within volume data range.")
|
||||
volume = np.array(volume, dtype=np.float64, order="C")
|
||||
|
||||
# Extract raw triangles using marching cubes in Cython
|
||||
# Returns a list of length-3 lists, each sub-list containing three
|
||||
# tuples. The tuples hold (x, y, z) coordinates for triangle vertices.
|
||||
# Note: this algorithm is fast, but returns degenerate "triangles" which
|
||||
# have repeated vertices - and equivalent vertices are redundantly
|
||||
# placed in every triangle they connect with.
|
||||
raw_tris = _marching_cubes_cy.iterate_and_store_3d(volume, float(level),
|
||||
sampling)
|
||||
|
||||
# Find and collect unique vertices, storing triangle verts as indices.
|
||||
# Returns a true mesh with no degenerate faces.
|
||||
verts, faces = _marching_cubes_cy.unpack_unique_verts(raw_tris)
|
||||
|
||||
return np.asarray(verts), np.asarray(faces)
|
||||
|
||||
|
||||
def mesh_surface_area(verts, tris):
|
||||
"""
|
||||
Compute surface area, given vertices & triangular faces
|
||||
|
||||
Parameters
|
||||
----------
|
||||
verts : (V, 3) array of floats
|
||||
Array containing (x, y, z) coordinates for V unique mesh vertices.
|
||||
faces : (F, 3) array of ints
|
||||
List of length-3 lists of integers, referencing vertex coordinates as
|
||||
provided in `verts`
|
||||
|
||||
Returns
|
||||
-------
|
||||
area : float
|
||||
Surface area of mesh. Units now [coordinate units] ** 2.
|
||||
|
||||
Notes
|
||||
-----
|
||||
The arguments expected by this function are the exact outputs from
|
||||
`skimage.measure.marching_cubes`. For unit correct output, ensure correct
|
||||
`spacing` was passed to `skimage.measure.marching_cubes`.
|
||||
|
||||
This algorithm works properly only if the ``faces`` provided are all
|
||||
triangles.
|
||||
|
||||
See Also
|
||||
--------
|
||||
skimage.measure.marching_cubes
|
||||
|
||||
"""
|
||||
# Fancy indexing to define two vector arrays from triangle vertices
|
||||
actual_verts = verts[tris]
|
||||
a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
|
||||
b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
|
||||
del actual_verts
|
||||
|
||||
# Area of triangle in 3D = 1/2 * Euclidean norm of cross product
|
||||
return ((np.cross(a, b) ** 2).sum(axis=1) ** 0.5).sum() / 2.
|
||||
@@ -0,0 +1,987 @@
|
||||
#cython: cdivision=True
|
||||
#cython: boundscheck=False
|
||||
#cython: nonecheck=False
|
||||
#cython: wraparound=False
|
||||
import numpy as np
|
||||
cimport numpy as cnp
|
||||
|
||||
|
||||
cdef inline double _get_fraction(double from_value, double to_value,
|
||||
double level):
|
||||
if (to_value == from_value):
|
||||
return 0
|
||||
return ((level - from_value) / (to_value - from_value))
|
||||
|
||||
|
||||
def unpack_unique_verts(list trilist):
|
||||
"""
|
||||
Convert a list of lists of tuples corresponding to triangle vertices
|
||||
into a unique vertex list, and a list of triangle faces w/indices
|
||||
corresponding to entries of the vertex list.
|
||||
|
||||
"""
|
||||
cdef Py_ssize_t idx = 0
|
||||
cdef Py_ssize_t n_tris = len(trilist)
|
||||
cdef Py_ssize_t i, j
|
||||
cdef dict vert_index = {}
|
||||
cdef list vert_list = []
|
||||
cdef list face_list = []
|
||||
cdef list templist
|
||||
|
||||
# Iterate over triangles
|
||||
for i in range(n_tris):
|
||||
templist = []
|
||||
|
||||
# Only parse vertices from non-degenerate triangles
|
||||
if not ((trilist[i][0] == trilist[i][1]) or
|
||||
(trilist[i][0] == trilist[i][2]) or
|
||||
(trilist[i][1] == trilist[i][2])):
|
||||
|
||||
# Iterate over vertices within each triangle
|
||||
for j in range(3):
|
||||
vert = trilist[i][j]
|
||||
|
||||
# Check if a new unique vertex found
|
||||
if vert not in vert_index:
|
||||
vert_index[vert] = idx
|
||||
templist.append(idx)
|
||||
vert_list.append(vert)
|
||||
idx += 1
|
||||
else:
|
||||
templist.append(vert_index[vert])
|
||||
|
||||
face_list.append(templist)
|
||||
|
||||
return vert_list, face_list
|
||||
|
||||
|
||||
def iterate_and_store_3d(double[:, :, ::1] arr, double level,
|
||||
tuple sampling=(1., 1., 1.)):
|
||||
"""Iterate across the given array in a marching-cubes fashion,
|
||||
looking for volumes with edges that cross 'level'. If such a volume is
|
||||
found, appropriate triangulations are added to a growing list of
|
||||
faces to be returned by this function.
|
||||
|
||||
If `sampling` is not provided, vertices are returned in the indexing
|
||||
coordinate system (assuming all 3 spatial dimensions sampled equally).
|
||||
If `sampling` is provided, vertices will be returned in volume coordinates
|
||||
relative to the origin, regularly spaced as specified in each dimension.
|
||||
|
||||
"""
|
||||
if arr.shape[0] < 2 or arr.shape[1] < 2 or arr.shape[2] < 2:
|
||||
raise ValueError("Input array must be at least 2x2x2.")
|
||||
if len(sampling) != 3:
|
||||
raise ValueError("`sampling` must be (double, double, double)")
|
||||
|
||||
cdef list face_list = []
|
||||
cdef list norm_list = []
|
||||
cdef Py_ssize_t n
|
||||
cdef bint odd_sampling, plus_z
|
||||
plus_z = False
|
||||
if [float(i) for i in sampling] == [1.0, 1.0, 1.0]:
|
||||
odd_sampling = False
|
||||
else:
|
||||
odd_sampling = True
|
||||
|
||||
# The plan is to iterate a 2x2x2 cube across the input array. This means
|
||||
# the upper-left corner of the cube needs to iterate across a sub-array
|
||||
# of size one-less-large in each direction (so we can get away with no
|
||||
# bounds checking in Cython). The cube is represented by eight vertices:
|
||||
# v1, v2, ..., v8, oriented thus (see Lorensen, Figure 4):
|
||||
#
|
||||
# v8 ------ v7
|
||||
# / | / | y
|
||||
# / | / | ^ z
|
||||
# v4 ------ v3 | | /
|
||||
# | v5 ----|- v6 |/ (note: NOT right handed!)
|
||||
# | / | / ----> x
|
||||
# | / | /
|
||||
# v1 ------ v2
|
||||
#
|
||||
# We also maintain the current 2D coordinates for v1, and ensure the array
|
||||
# is of type 'double' and is C-contiguous (last index varies fastest).
|
||||
|
||||
# Coords start at (0, 0, 0).
|
||||
cdef Py_ssize_t[3] coords
|
||||
coords[0] = 0
|
||||
coords[1] = 0
|
||||
coords[2] = 0
|
||||
|
||||
# Extract doubles from `sampling` for speed
|
||||
cdef double[3] sampling2
|
||||
sampling2[0] = sampling[0]
|
||||
sampling2[1] = sampling[1]
|
||||
sampling2[2] = sampling[2]
|
||||
|
||||
# Calculate the number of iterations we'll need
|
||||
cdef Py_ssize_t num_cube_steps = ((arr.shape[0] - 1) *
|
||||
(arr.shape[1] - 1) *
|
||||
(arr.shape[2] - 1))
|
||||
|
||||
cdef unsigned char cube_case = 0
|
||||
cdef tuple e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12
|
||||
cdef double v1, v2, v3, v4, v5, v6, v7, v8, r0, r1, c0, c1, d0, d1
|
||||
cdef Py_ssize_t x0, y0, z0, x1, y1, z1
|
||||
e5, e6, e7, e8 = (0, 0, 0), (0, 0, 0), (0, 0, 0), (0, 0, 0)
|
||||
|
||||
for n in range(num_cube_steps):
|
||||
# There are 255 unique values for `cube_case`. This algorithm follows
|
||||
# the Lorensen paper in vertex and edge labeling, however, it should
|
||||
# be noted that Lorensen used a left-handed coordinate system while
|
||||
# NumPy uses a proper right handed system. Transforming between these
|
||||
# coordinate systems was handled in the definitions of the cube
|
||||
# vertices v1, v2, ..., v8.
|
||||
#
|
||||
# Refer to the paper, figure 4, for cube edge designations e1, ... e12
|
||||
|
||||
# Standard Py_ssize_t coordinates for indexing
|
||||
x0, y0, z0 = coords[0], coords[1], coords[2]
|
||||
x1, y1, z1 = x0 + 1, y0 + 1, z0 + 1
|
||||
|
||||
if odd_sampling:
|
||||
# These doubles are the modified world coordinates; they are only
|
||||
# calculated if non-default `sampling` provided.
|
||||
r0 = coords[0] * sampling2[0]
|
||||
c0 = coords[1] * sampling2[1]
|
||||
d0 = coords[2] * sampling2[2]
|
||||
r1 = r0 + sampling2[0]
|
||||
c1 = c0 + sampling2[1]
|
||||
d1 = d0 + sampling2[2]
|
||||
else:
|
||||
r0, c0, d0, r1, c1, d1 = x0, y0, z0, x1, y1, z1
|
||||
|
||||
# We use a right-handed coordinate system, UNlike the paper, but want
|
||||
# to index in agreement - the coordinate adjustment takes place here.
|
||||
v1 = arr[x0, y0, z0]
|
||||
v2 = arr[x1, y0, z0]
|
||||
v3 = arr[x1, y1, z0]
|
||||
v4 = arr[x0, y1, z0]
|
||||
v5 = arr[x0, y0, z1]
|
||||
v6 = arr[x1, y0, z1]
|
||||
v7 = arr[x1, y1, z1]
|
||||
v8 = arr[x0, y1, z1]
|
||||
|
||||
# Unique triangulation cases
|
||||
cube_case = 0
|
||||
if (v1 > level): cube_case += 1
|
||||
if (v2 > level): cube_case += 2
|
||||
if (v3 > level): cube_case += 4
|
||||
if (v4 > level): cube_case += 8
|
||||
if (v5 > level): cube_case += 16
|
||||
if (v6 > level): cube_case += 32
|
||||
if (v7 > level): cube_case += 64
|
||||
if (v8 > level): cube_case += 128
|
||||
|
||||
if (cube_case != 0 and cube_case != 255):
|
||||
# Only do anything if there's a plane intersecting the cube.
|
||||
# Cases 0 and 255 are entirely below/above the contour.
|
||||
|
||||
if cube_case > 127:
|
||||
if ((cube_case != 150) and
|
||||
(cube_case != 170) and
|
||||
(cube_case != 195)):
|
||||
cube_case = 255 - cube_case
|
||||
|
||||
# Calculate cube edges, to become triangulation vertices.
|
||||
# If we moved in a convenient direction, save 1/3 of the effort by
|
||||
# re-assigning prior results.
|
||||
if plus_z:
|
||||
# Reassign prior calculated edges
|
||||
e1 = e5
|
||||
e2 = e6
|
||||
e3 = e7
|
||||
e4 = e8
|
||||
else:
|
||||
# Calculate edges normally
|
||||
if odd_sampling:
|
||||
e1 = r0 + _get_fraction(v1, v2, level) * sampling2[0], c0, d0
|
||||
e2 = r1, c0 + _get_fraction(v2, v3, level) * sampling2[1], d0
|
||||
e3 = r0 + _get_fraction(v4, v3, level) * sampling2[0], c1, d0
|
||||
e4 = r0, c0 + _get_fraction(v1, v4, level) * sampling2[1], d0
|
||||
else:
|
||||
e1 = r0 + _get_fraction(v1, v2, level), c0, d0
|
||||
e2 = r1, c0 + _get_fraction(v2, v3, level), d0
|
||||
e3 = r0 + _get_fraction(v4, v3, level), c1, d0
|
||||
e4 = r0, c0 + _get_fraction(v1, v4, level), d0
|
||||
|
||||
# These must be calculated at each point unless we implemented a
|
||||
# large, growing lookup table for all adjacent values; could save
|
||||
# ~30% in terms of runtime at the expense of memory usage and
|
||||
# much greater complexity.
|
||||
if odd_sampling:
|
||||
e5 = r0 + _get_fraction(v5, v6, level) * sampling2[0], c0, d1
|
||||
e6 = r1, c0 + _get_fraction(v6, v7, level) * sampling2[1], d1
|
||||
e7 = r0 + _get_fraction(v8, v7, level) * sampling2[0], c1, d1
|
||||
e8 = r0, c0 + _get_fraction(v5, v8, level) * sampling2[1], d1
|
||||
e9 = r0, c0, d0 + _get_fraction(v1, v5, level) * sampling2[2]
|
||||
e10 = r1, c0, d0 + _get_fraction(v2, v6, level) * sampling2[2]
|
||||
e11 = r0, c1, d0 + _get_fraction(v4, v8, level) * sampling2[2]
|
||||
e12 = r1, c1, d0 + _get_fraction(v3, v7, level) * sampling2[2]
|
||||
else:
|
||||
e5 = r0 + _get_fraction(v5, v6, level), c0, d1
|
||||
e6 = r1, c0 + _get_fraction(v6, v7, level), d1
|
||||
e7 = r0 + _get_fraction(v8, v7, level), c1, d1
|
||||
e8 = r0, c0 + _get_fraction(v5, v8, level), d1
|
||||
e9 = r0, c0, d0 + _get_fraction(v1, v5, level)
|
||||
e10 = r1, c0, d0 + _get_fraction(v2, v6, level)
|
||||
e11 = r0, c1, d0 + _get_fraction(v4, v8, level)
|
||||
e12 = r1, c1, d0 + _get_fraction(v3, v7, level)
|
||||
|
||||
|
||||
# Append appropriate triangles to the growing output `face_list`
|
||||
_append_tris(face_list, cube_case, e1, e2, e3, e4, e5,
|
||||
e6, e7, e8, e9, e10, e11, e12)
|
||||
|
||||
# Advance the coords indices
|
||||
if coords[2] < arr.shape[2] - 2:
|
||||
coords[2] += 1
|
||||
plus_z = True
|
||||
elif coords[1] < arr.shape[1] - 2:
|
||||
coords[1] += 1
|
||||
coords[2] = 0
|
||||
plus_z = False
|
||||
else:
|
||||
coords[0] += 1
|
||||
coords[1] = 0
|
||||
coords[2] = 0
|
||||
plus_z = False
|
||||
|
||||
return face_list
|
||||
|
||||
|
||||
def _append_tris(list face_list, unsigned char case, tuple e1, tuple e2,
|
||||
tuple e3, tuple e4, tuple e5, tuple e6, tuple e7, tuple e8,
|
||||
tuple e9, tuple e10, tuple e11, tuple e12):
|
||||
# Permits recursive use for duplicated planes to conserve code - it's
|
||||
# quite long enough as-is.
|
||||
|
||||
if (case == 1):
|
||||
# front lower left corner
|
||||
face_list.append([e1, e4, e9])
|
||||
elif (case == 2):
|
||||
# front lower right corner
|
||||
face_list.append([e10, e2, e1])
|
||||
elif (case == 3):
|
||||
# front lower plane
|
||||
face_list.append([e2, e4, e9])
|
||||
face_list.append([e2, e9, e10])
|
||||
elif (case == 4):
|
||||
# front upper right corner
|
||||
face_list.append([e12, e3, e2])
|
||||
elif (case == 5):
|
||||
# lower left, upper right corners
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 6):
|
||||
# front right plane
|
||||
face_list.append([e12, e3, e1])
|
||||
face_list.append([e12, e1, e10])
|
||||
elif (case == 7):
|
||||
# Shelf including v1, v2, v3
|
||||
face_list.append([e3, e4, e12])
|
||||
face_list.append([e4, e9, e12])
|
||||
face_list.append([e12, e9, e10])
|
||||
elif (case == 8):
|
||||
# front upper left corner
|
||||
face_list.append([e3, e11, e4])
|
||||
elif (case == 9):
|
||||
# front left plane
|
||||
face_list.append([e3, e11, e9])
|
||||
face_list.append([e3, e9, e1])
|
||||
elif (case == 10):
|
||||
# upper left, lower right corners
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 11):
|
||||
# Shelf including v4, v1, v2
|
||||
face_list.append([e3, e11, e2])
|
||||
face_list.append([e11, e10, e2])
|
||||
face_list.append([e11, e9, e10])
|
||||
elif (case == 12):
|
||||
# front upper plane
|
||||
face_list.append([e11, e4, e12])
|
||||
face_list.append([e2, e4, e12])
|
||||
elif (case == 13):
|
||||
# Shelf including v1, v4, v3
|
||||
face_list.append([e11, e9, e12])
|
||||
face_list.append([e12, e9, e1])
|
||||
face_list.append([e12, e1, e2])
|
||||
elif (case == 14):
|
||||
# Shelf including v2, v3, v4
|
||||
face_list.append([e11, e10, e12])
|
||||
face_list.append([e11, e4, e10])
|
||||
face_list.append([e4, e1, e10])
|
||||
elif (case == 15):
|
||||
# Plane parallel to x-axis through middle
|
||||
face_list.append([e11, e9, e12])
|
||||
face_list.append([e12, e9, e10])
|
||||
elif (case == 16):
|
||||
# back lower left corner
|
||||
face_list.append([e8, e9, e5])
|
||||
elif (case == 17):
|
||||
# lower left plane
|
||||
face_list.append([e4, e1, e8])
|
||||
face_list.append([e8, e1, e5])
|
||||
elif (case == 18):
|
||||
# lower left back, lower right front corners
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 19):
|
||||
# Shelf including v1, v2, v5
|
||||
face_list.append([e8, e4, e2])
|
||||
face_list.append([e8, e2, e10])
|
||||
face_list.append([e8, e10, e5])
|
||||
elif (case == 20):
|
||||
# lower left back, upper right front corners
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 21):
|
||||
# lower left plane + upper right front corner, v1, v3, v5
|
||||
_append_tris(face_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 22):
|
||||
# front right plane + lower left back corner, v2, v3, v5
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 23):
|
||||
# Rotated case 14 in the paper
|
||||
face_list.append([e3, e10, e8])
|
||||
face_list.append([e3, e10, e12])
|
||||
face_list.append([e8, e10, e5])
|
||||
face_list.append([e3, e4, e8])
|
||||
elif (case == 24):
|
||||
# upper front left, lower back left corners
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 25):
|
||||
# Shelf including v1, v4, v5
|
||||
face_list.append([e1, e5, e3])
|
||||
face_list.append([e3, e8, e11])
|
||||
face_list.append([e3, e5, e8])
|
||||
elif (case == 26):
|
||||
# Three isolated corners
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 27):
|
||||
# Full corner v1, case 9 in paper: (v1, v2, v4, v5)
|
||||
face_list.append([e11, e3, e2])
|
||||
face_list.append([e11, e2, e10])
|
||||
face_list.append([e10, e11, e8])
|
||||
face_list.append([e8, e5, e10])
|
||||
elif (case == 28):
|
||||
# upper front plane + corner v5
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 29):
|
||||
# special case of 11 in the paper: (v1, v3, v4, v5)
|
||||
face_list.append([e11, e5, e2])
|
||||
face_list.append([e11, e12, e2])
|
||||
face_list.append([e11, e5, e8])
|
||||
face_list.append([e2, e1, e5])
|
||||
elif (case == 30):
|
||||
# Shelf (v2, v3, v4) and lower left back corner
|
||||
_append_tris(face_list, 14, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 31):
|
||||
# Shelf: (v6, v7, v8) by inversion
|
||||
face_list.append([e11, e12, e10])
|
||||
face_list.append([e11, e8, e10])
|
||||
face_list.append([e8, e10, e5])
|
||||
elif (case == 32):
|
||||
# lower right back corner
|
||||
face_list.append([e6, e5, e10])
|
||||
elif (case == 33):
|
||||
# lower right back, lower left front corners
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 34):
|
||||
# lower right plane
|
||||
face_list.append([e1, e2, e5])
|
||||
face_list.append([e2, e6, e5])
|
||||
elif (case == 35):
|
||||
# Shelf: v1, v2, v6
|
||||
face_list.append([e4, e2, e6])
|
||||
face_list.append([e4, e9, e6])
|
||||
face_list.append([e6, e9, e5])
|
||||
elif (case == 36):
|
||||
# upper right front, lower right back corners
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 37):
|
||||
# lower left front, upper right front, lower right back corners
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 38):
|
||||
# Shelf: v2, v3, v6
|
||||
face_list.append([e3, e1, e5])
|
||||
face_list.append([e3, e5, e12])
|
||||
face_list.append([e12, e5, e6])
|
||||
elif (case == 39):
|
||||
# Full corner v2: (v1, v2, v3, v6)
|
||||
face_list.append([e3, e4, e5])
|
||||
face_list.append([e4, e9, e5])
|
||||
face_list.append([e3, e5, e6])
|
||||
face_list.append([e3, e12, e6])
|
||||
elif (case == 40):
|
||||
# upper left front, lower right back corners
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 41):
|
||||
# front left plane, lower right back corner
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 42):
|
||||
# lower right plane, upper front left corner
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 43):
|
||||
# Rotated case 11 in paper
|
||||
face_list.append([e11, e3, e9])
|
||||
face_list.append([e3, e9, e6])
|
||||
face_list.append([e3, e2, e6])
|
||||
face_list.append([e9, e5, e6])
|
||||
elif (case == 44):
|
||||
# upper front plane, lower right back corner
|
||||
_append_tris(face_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 45):
|
||||
# Shelf: (v1, v3, v4) + lower right back corner
|
||||
_append_tris(face_list, 13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 46):
|
||||
# Rotated case 14 in paper
|
||||
face_list.append([e4, e11, e12])
|
||||
face_list.append([e4, e12, e5])
|
||||
face_list.append([e12, e5, e6])
|
||||
face_list.append([e4, e5, e1])
|
||||
elif (case == 47):
|
||||
# Shelf: (v5, v8, v7) by inversion
|
||||
face_list.append([e11, e9, e12])
|
||||
face_list.append([e12, e9, e5])
|
||||
face_list.append([e12, e5, e6])
|
||||
elif (case == 48):
|
||||
# Back lower plane
|
||||
face_list.append([e9, e10, e6])
|
||||
face_list.append([e9, e6, e8])
|
||||
elif (case == 49):
|
||||
# Shelf: (v1, v5, v6)
|
||||
face_list.append([e4, e8, e6])
|
||||
face_list.append([e4, e6, e1])
|
||||
face_list.append([e6, e1, e10])
|
||||
elif (case == 50):
|
||||
# Shelf: (v2, v5, v6)
|
||||
face_list.append([e8, e6, e2])
|
||||
face_list.append([e8, e2, e1])
|
||||
face_list.append([e8, e9, e1])
|
||||
elif (case == 51):
|
||||
# Plane through middle of cube, parallel to x-z axis
|
||||
face_list.append([e4, e8, e2])
|
||||
face_list.append([e8, e2, e6])
|
||||
elif (case == 52):
|
||||
# Back lower plane, and front upper right corner
|
||||
_append_tris(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 53):
|
||||
# Shelf (v1, v5, v6) and front upper right corner
|
||||
_append_tris(face_list, 49, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 54):
|
||||
# Rotated case 11 from paper (v2, v3, v5, v6)
|
||||
face_list.append([e1, e9, e3])
|
||||
face_list.append([e9, e3, e6])
|
||||
face_list.append([e9, e8, e6])
|
||||
face_list.append([e12, e3, e6])
|
||||
elif (case == 55):
|
||||
# Shelf: (v4, v8, v7) by inversion
|
||||
face_list.append([e4, e8, e6])
|
||||
face_list.append([e4, e6, e3])
|
||||
face_list.append([e6, e3, e12])
|
||||
elif (case == 56):
|
||||
# Back lower plane + upper left front corner
|
||||
_append_tris(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 57):
|
||||
# Rotated case 14 from paper (v4, v1, v5, v6)
|
||||
face_list.append([e3, e11, e8])
|
||||
face_list.append([e3, e8, e10])
|
||||
face_list.append([e10, e6, e8])
|
||||
face_list.append([e3, e1, e10])
|
||||
elif (case == 58):
|
||||
# Shelf: (v2, v6, v5) + upper left front corner
|
||||
_append_tris(face_list, 50, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 59):
|
||||
# Shelf: (v3, v7, v8) by inversion
|
||||
face_list.append([e2, e6, e8])
|
||||
face_list.append([e8, e2, e3])
|
||||
face_list.append([e8, e3, e11])
|
||||
elif (case == 60):
|
||||
# AMBIGUOUS CASE: parallel planes (front upper, back lower)
|
||||
_append_tris(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 61):
|
||||
# Upper back plane + lower right front corner by inversion
|
||||
_append_tris(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 62):
|
||||
# Upper back plane + lower left front corner by inversion
|
||||
_append_tris(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 63):
|
||||
# Upper back plane
|
||||
face_list.append([e11, e12, e6])
|
||||
face_list.append([e11, e8, e6])
|
||||
elif (case == 64):
|
||||
# Upper right back corner
|
||||
face_list.append([e12, e7, e6])
|
||||
elif (case == 65):
|
||||
# upper right back, lower left front corners
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 66):
|
||||
# upper right back, lower right front corners
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 67):
|
||||
# lower front plane + upper right back corner
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 68):
|
||||
# upper right plane
|
||||
face_list.append([e3, e2, e6])
|
||||
face_list.append([e3, e7, e6])
|
||||
elif (case == 69):
|
||||
# Upper right plane, lower left front corner
|
||||
_append_tris(face_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 70):
|
||||
# Shelf: (v2, v3, v7)
|
||||
face_list.append([e1, e3, e7])
|
||||
face_list.append([e1, e10, e7])
|
||||
face_list.append([e7, e10, e6])
|
||||
elif (case == 71):
|
||||
# Rotated version of case 11 in paper (v1, v2, v3, v7)
|
||||
face_list.append([e10, e7, e4])
|
||||
face_list.append([e4, e3, e7])
|
||||
face_list.append([e10, e4, e9])
|
||||
face_list.append([e7, e10, e6])
|
||||
elif (case == 72):
|
||||
# upper left front, upper right back corners
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 73):
|
||||
# front left plane, upper right back corner
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 74):
|
||||
# Three isolated corners, exactly case 7 in paper
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 75):
|
||||
# Shelf: (v1, v2, v4) + upper right back corner
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 11, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 76):
|
||||
# Shelf: (v4, v3, v7)
|
||||
face_list.append([e4, e2, e6])
|
||||
face_list.append([e4, e11, e7])
|
||||
face_list.append([e4, e7, e6])
|
||||
elif (case == 77):
|
||||
# Rotated case 14 in paper (v1, v4, v3, v7)
|
||||
face_list.append([e11, e9, e1])
|
||||
face_list.append([e11, e1, e6])
|
||||
face_list.append([e1, e6, e2])
|
||||
face_list.append([e11, e6, e7])
|
||||
elif (case == 78):
|
||||
# Full corner v3: (v2, v3, v4, v7)
|
||||
face_list.append([e1, e4, e7])
|
||||
face_list.append([e1, e7, e6])
|
||||
face_list.append([e4, e11, e7])
|
||||
face_list.append([e1, e10, e6])
|
||||
elif (case == 79):
|
||||
# Shelf: (v6, v5, v8) by inversion
|
||||
face_list.append([e9, e11, e10])
|
||||
face_list.append([e11, e7, e10])
|
||||
face_list.append([e7, e10, e6])
|
||||
elif (case == 80):
|
||||
# lower left back, upper right back corners (v5, v7)
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 81):
|
||||
# lower left plane, upper right back corner
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 82):
|
||||
# isolated corners (v2, v5, v7)
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 83):
|
||||
# Shelf: (v1, v2, v5) + upper right back corner
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 19, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 84):
|
||||
# upper right plane, lower left back corner
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 85):
|
||||
# AMBIGUOUS CASE: upper right and lower left parallel planes
|
||||
_append_tris(face_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 86):
|
||||
# Shelf: (v2, v3, v7) + lower left back corner
|
||||
_append_tris(face_list, 70, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 87):
|
||||
# Upper left plane + lower right back corner, by inversion
|
||||
_append_tris(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 88):
|
||||
# Isolated corners v4, v5, v7
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 89):
|
||||
# Shelf: (v1, v4, v5) + isolated corner v7
|
||||
_append_tris(face_list, 25, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 90):
|
||||
# Four isolated corners v2, v4, v5, v7
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 91):
|
||||
# Three isolated corners, v3, v6, v8 by inversion
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 92):
|
||||
# Shelf (v4, v3, v7) + isolated corner v5
|
||||
_append_tris(face_list, 76, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 93):
|
||||
# Lower right plane + isolated corner v8 by inversion
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 94):
|
||||
# Isolated corners v1, v6, v8 by inversion
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 95):
|
||||
# Isolated corners v6, v8 by inversion
|
||||
_append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 96):
|
||||
# back right plane
|
||||
face_list.append([e7, e12, e5])
|
||||
face_list.append([e5, e10, e12])
|
||||
elif (case == 97):
|
||||
# back right plane + isolated corner v1
|
||||
_append_tris(face_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 98):
|
||||
# Shelf: (v2, v6, v7)
|
||||
face_list.append([e1, e7, e5])
|
||||
face_list.append([e7, e1, e12])
|
||||
face_list.append([e1, e12, e2])
|
||||
elif (case == 99):
|
||||
# Rotated case 14 in paper: (v1, v2, v6, v7)
|
||||
face_list.append([e9, e2, e7])
|
||||
face_list.append([e9, e2, e4])
|
||||
face_list.append([e2, e7, e12])
|
||||
face_list.append([e7, e9, e5])
|
||||
elif (case == 100):
|
||||
# Shelf: (v3, v6, v7)
|
||||
face_list.append([e3, e7, e5])
|
||||
face_list.append([e3, e5, e2])
|
||||
face_list.append([e2, e5, e10])
|
||||
elif (case == 101):
|
||||
# Shelf: (v3, v6, v7) + isolated corner v1
|
||||
_append_tris(face_list, 100, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 102):
|
||||
# Plane bisecting left-right halves of cube
|
||||
face_list.append([e1, e3, e7])
|
||||
face_list.append([e1, e7, e5])
|
||||
elif (case == 103):
|
||||
# Shelf: (v4, v5, v8) by inversion
|
||||
face_list.append([e3, e7, e5])
|
||||
face_list.append([e3, e5, e4])
|
||||
face_list.append([e4, e5, e9])
|
||||
elif (case == 104):
|
||||
# Back right plane + isolated corner v4
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 105):
|
||||
# AMBIGUOUS CASE: back right and front left planes
|
||||
_append_tris(face_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 106):
|
||||
# Shelf: (v2, v6, v7) + isolated corner v4
|
||||
_append_tris(face_list, 98, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 107):
|
||||
# Back left plane + isolated corner v3 by inversion
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 108):
|
||||
# Rotated case 11 from paper: (v4, v3, v7, v6)
|
||||
face_list.append([e4, e10, e7])
|
||||
face_list.append([e4, e10, e2])
|
||||
face_list.append([e4, e11, e7])
|
||||
face_list.append([e7, e10, e5])
|
||||
elif (case == 109):
|
||||
# Back left plane + isolated corner v2 by inversion
|
||||
_append_tris(face_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 110):
|
||||
# Shelf: (v1, v5, v8) by inversion
|
||||
face_list.append([e1, e5, e7])
|
||||
face_list.append([e1, e7, e11])
|
||||
face_list.append([e1, e11, e4])
|
||||
elif (case == 111):
|
||||
# Back left plane
|
||||
face_list.append([e11, e9, e7])
|
||||
face_list.append([e9, e7, e5])
|
||||
elif (case == 112):
|
||||
# Shelf: (v5, v6, v7)
|
||||
face_list.append([e9, e10, e12])
|
||||
face_list.append([e9, e12, e7])
|
||||
face_list.append([e9, e7, e8])
|
||||
elif (case == 113):
|
||||
# Exactly case 11 from paper: (v1, v5, v6, v7)
|
||||
face_list.append([e1, e8, e12])
|
||||
face_list.append([e1, e8, e4])
|
||||
face_list.append([e8, e7, e12])
|
||||
face_list.append([e12, e1, e10])
|
||||
elif (case == 114):
|
||||
# Full corner v6: (v2, v6, v7, v5)
|
||||
face_list.append([e1, e9, e7])
|
||||
face_list.append([e1, e7, e12])
|
||||
face_list.append([e1, e12, e2])
|
||||
face_list.append([e9, e8, e7])
|
||||
elif (case == 115):
|
||||
# Shelf: (v3, v4, v8)
|
||||
face_list.append([e2, e4, e8])
|
||||
face_list.append([e2, e12, e7])
|
||||
face_list.append([e2, e8, e7])
|
||||
elif (case == 116):
|
||||
# Rotated case 14 in paper: (v5, v6, v7, v3)
|
||||
face_list.append([e9, e2, e7])
|
||||
face_list.append([e9, e2, e10])
|
||||
face_list.append([e9, e8, e7])
|
||||
face_list.append([e2, e3, e7])
|
||||
elif (case == 117):
|
||||
# upper left plane + isolated corner v2 by inversion
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 118):
|
||||
# Shelf: (v1, v4, v8)
|
||||
face_list.append([e1, e3, e7])
|
||||
face_list.append([e7, e1, e8])
|
||||
face_list.append([e1, e8, e9])
|
||||
elif (case == 119):
|
||||
# Upper left plane
|
||||
face_list.append([e4, e3, e7])
|
||||
face_list.append([e4, e8, e7])
|
||||
elif (case == 120):
|
||||
# Shelf: (v5, v6, v7) + isolated corner v4
|
||||
_append_tris(face_list, 112, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 121):
|
||||
# Front right plane + isolated corner v8
|
||||
_append_tris(face_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 122):
|
||||
# Isolated corners v1, v3, v8
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 123):
|
||||
# Isolated corners v3, v8
|
||||
_append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 124):
|
||||
# Front lower plane + isolated corner v8
|
||||
_append_tris(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 125):
|
||||
# Isolated corners v2, v8
|
||||
_append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 126):
|
||||
# Isolated corners v1, v8
|
||||
_append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 127):
|
||||
# Isolated corner v8
|
||||
face_list.append([e11, e7, e8])
|
||||
elif (case == 150):
|
||||
# AMBIGUOUS CASE: back right and front left planes
|
||||
# In these cube_case > 127 cases, the vertices are identical BUT
|
||||
# they are connected in the opposite fashion.
|
||||
_append_tris(face_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 170):
|
||||
# AMBIGUOUS CASE: upper left and lower right planes
|
||||
# In these cube_case > 127 cases, the vertices are identical BUT
|
||||
# they are connected in the opposite fashion.
|
||||
_append_tris(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
elif (case == 195):
|
||||
# AMBIGUOUS CASE: back upper and front lower planes
|
||||
# In these cube_case > 127 cases, the vertices are identical BUT
|
||||
# they are connected in the opposite fashion.
|
||||
_append_tris(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
_append_tris(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
|
||||
e11, e12)
|
||||
|
||||
return
|
||||
@@ -14,11 +14,15 @@ def configuration(parent_package='', top_path=None):
|
||||
|
||||
cython(['_find_contours.pyx'], working_path=base_path)
|
||||
cython(['_moments.pyx'], working_path=base_path)
|
||||
cython(['_marching_cubes_cy.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_find_contours', sources=['_find_contours.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_moments', sources=['_moments.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_marching_cubes_cy',
|
||||
sources=['_marching_cubes_cy.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -0,0 +1,40 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_raises
|
||||
|
||||
from skimage.draw import ellipsoid, ellipsoid_stats
|
||||
from skimage.measure import marching_cubes, mesh_surface_area
|
||||
|
||||
|
||||
def test_marching_cubes_isotropic():
|
||||
ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True)
|
||||
_, surf = ellipsoid_stats(6, 10, 16)
|
||||
verts, faces = marching_cubes(ellipsoid_isotropic, 0.)
|
||||
surf_calc = mesh_surface_area(verts, faces)
|
||||
|
||||
# Test within 1% tolerance for isotropic. Will always underestimate.
|
||||
assert surf > surf_calc and surf_calc > surf * 0.99
|
||||
|
||||
|
||||
def test_marching_cubes_anisotropic():
|
||||
sampling = (1., 10 / 6., 16 / 6.)
|
||||
ellipsoid_anisotropic = ellipsoid(6, 10, 16, sampling=sampling,
|
||||
levelset=True)
|
||||
_, surf = ellipsoid_stats(6, 10, 16, sampling=sampling)
|
||||
verts, faces = marching_cubes(ellipsoid_anisotropic, 0.,
|
||||
sampling=sampling)
|
||||
surf_calc = mesh_surface_area(verts, faces)
|
||||
|
||||
# Test within 1.5% tolerance for anisotropic. Will always underestimate.
|
||||
assert surf > surf_calc and surf_calc > surf * 0.985
|
||||
|
||||
|
||||
def test_invalid_input():
|
||||
assert_raises(ValueError, marching_cubes, np.zeros((2, 2, 1)), 0)
|
||||
assert_raises(ValueError, marching_cubes, np.zeros((2, 2, 1)), 1)
|
||||
assert_raises(ValueError, marching_cubes, np.ones((3, 3, 3)), 1,
|
||||
sampling=(1, 2))
|
||||
assert_raises(ValueError, marching_cubes, np.zeros((20, 20)), 0)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
np.testing.run_module_suite()
|
||||
Reference in New Issue
Block a user