MarchingCases
vtk-examples/Python/VisualizationAlgorithms/MarchingCases
Description¶
This example will help you understand the Marching Cubes Algorithm. The example takes one optional argument, a case number. There are 15 Marching Cubes cases, 0-14. There are also 15 complementary cases where the inside/outside value is flipped. To see a complementary case, supply a negative case number. For example, -7 is the complementary case of 7.
Note
According to the ACM Digital Library, the Marching Cubes paper is the most cited Siggraph paper.
Other languages
See (Cxx)
Question
If you have a question about this example, please use the VTK Discourse Forum
Code¶
MarchingCases.py
#!/usr/bin/env python
# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonCore import (
vtkFloatArray,
vtkIdList,
vtkPoints
)
from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid
from vtkmodules.vtkCommonTransforms import vtkTransform
from vtkmodules.vtkFiltersCore import (
vtkContourFilter,
vtkGlyph3D,
vtkThresholdPoints,
vtkTubeFilter
)
# vtkExtractEdges moved from vtkFiltersExtraction to vtkFiltersCore in
# VTK commit d9981b9aeb93b42d1371c6e295d76bfdc18430bd
try:
from vtkmodules.vtkFiltersCore import vtkExtractEdges
except ImportError:
from vtkmodules.vtkFiltersExtraction import vtkExtractEdges
from vtkmodules.vtkFiltersGeneral import (
vtkShrinkPolyData,
vtkTransformPolyDataFilter
)
from vtkmodules.vtkFiltersSources import (
vtkCubeSource,
vtkSphereSource
)
from vtkmodules.vtkRenderingCore import (
vtkActor,
vtkPolyDataMapper,
vtkRenderWindow,
vtkRenderWindowInteractor,
vtkRenderer
)
from vtkmodules.vtkRenderingFreeType import vtkVectorText
def main():
mc_cases, rotation, label = get_program_parameters()
if not mc_cases:
mc_cases = [7]
else:
# Ensure that they are unique.
mc_cases = list(set(mc_cases))
# Check that they lie in the correct range.
badCases = []
for item in mc_cases:
if abs(int(item) > 14):
badCases.append(item)
if badCases:
print('Bad case number(s)', ','.join(map(str, badCases)))
for item in badCases:
mc_cases.remove(item)
if not mc_cases:
print('No cases.')
return
marching_cubes(mc_cases, rotation, label)
def get_program_parameters():
import argparse
description = 'Marching cubes cases for 3D isosurface generation.'
epilogue = '''
Marching cubes cases for 3D isosurface generation.
The 256 possible cases have been reduced to 15 cases using symmetry.
Dark vertices are greater than the selected isosurface value.
For the cases, enter them as integers separated by a space e.g: 1 2 3
'''
parser = argparse.ArgumentParser(description=description, epilog=epilogue,
formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('cases', nargs='*', type=int, default=[],
help='A list of integers i such that 0 <= abs(i) < 14, corresponding to the cases desired.')
parser.add_argument('-r', '--rotation', type=int, default=0,
help='Rotate camera around the cube, for i such that 0 <= abs(i) < 4,\
corresponding to 0, 90, 180, 270 degrees.')
# Use a mutually exclusive group.
label_parser = parser.add_mutually_exclusive_group(required=False)
label_parser.add_argument('-l', '--label', action='store_true', dest='label',
help='Display a label, true by default.')
label_parser.add_argument('-n', '--no_label', action='store_false', dest='label',
help='Supress diaplaying a label.')
parser.set_defaults(label=True)
args = parser.parse_args()
return args.cases, args.rotation, args.label
def marching_cubes(mcCases, rotation=0, label=True):
color = vtkNamedColors()
# Rotate the final figure 0, 90, 180, 270 degrees.
rotation = abs(int(rotation))
if rotation > 3:
rotation = 0
if len(mcCases) > 1:
print('Cases', ', '.join(map(str, mcCases)))
else:
print('Cases', ','.join(map(str, mcCases)))
print('Rotated', rotation * 90, 'degrees.')
renWin = vtkRenderWindow()
renWin.SetSize(640, 480)
iren = vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
# Always use a grid of four columns unless number of cases < 4.
renderers = list()
gridSize = ((len(mcCases) + 3) // 4) * 4
if len(mcCases) < 4:
gridSize = len(mcCases)
for i in range(0, gridSize):
# Create the Renderer
renderer = vtkRenderer()
renderers.append(renderer)
# Set the background color.
renderers[i].SetBackground(color.GetColor3d('slate_grey'))
renWin.AddRenderer(renderer)
for i in range(0, len(mcCases)):
# Define a Single Cube
Scalars = vtkFloatArray()
Scalars.InsertNextValue(1.0)
Scalars.InsertNextValue(0.0)
Scalars.InsertNextValue(0.0)
Scalars.InsertNextValue(1.0)
Scalars.InsertNextValue(0.0)
Scalars.InsertNextValue(0.0)
Scalars.InsertNextValue(0.0)
Scalars.InsertNextValue(0.0)
Points = vtkPoints()
Points.InsertNextPoint(0, 0, 0)
Points.InsertNextPoint(1, 0, 0)
Points.InsertNextPoint(1, 1, 0)
Points.InsertNextPoint(0, 1, 0)
Points.InsertNextPoint(0, 0, 1)
Points.InsertNextPoint(1, 0, 1)
Points.InsertNextPoint(1, 1, 1)
Points.InsertNextPoint(0, 1, 1)
Ids = vtkIdList()
Ids.InsertNextId(0)
Ids.InsertNextId(1)
Ids.InsertNextId(2)
Ids.InsertNextId(3)
Ids.InsertNextId(4)
Ids.InsertNextId(5)
Ids.InsertNextId(6)
Ids.InsertNextId(7)
Grid = vtkUnstructuredGrid()
Grid.Allocate(10, 10)
Grid.InsertNextCell(12, Ids)
Grid.SetPoints(Points)
Grid.GetPointData().SetScalars(Scalars)
# Find the triangles that lie along the 0.5 contour in this cube.
Marching = vtkContourFilter()
Marching.SetInputData(Grid)
Marching.SetValue(0, 0.5)
Marching.Update()
# Extract the edges of the triangles just found.
triangleEdges = vtkExtractEdges()
triangleEdges.SetInputConnection(Marching.GetOutputPort())
# Draw the edges as tubes instead of lines. Also create the associated
# mapper and actor to display the tubes.
triangleEdgeTubes = vtkTubeFilter()
triangleEdgeTubes.SetInputConnection(triangleEdges.GetOutputPort())
triangleEdgeTubes.SetRadius(.005)
triangleEdgeTubes.SetNumberOfSides(6)
triangleEdgeTubes.UseDefaultNormalOn()
triangleEdgeTubes.SetDefaultNormal(.577, .577, .577)
triangleEdgeMapper = vtkPolyDataMapper()
triangleEdgeMapper.SetInputConnection(triangleEdgeTubes.GetOutputPort())
triangleEdgeMapper.ScalarVisibilityOff()
triangleEdgeActor = vtkActor()
triangleEdgeActor.SetMapper(triangleEdgeMapper)
triangleEdgeActor.GetProperty().SetDiffuseColor(
color.GetColor3d('lamp_black'))
triangleEdgeActor.GetProperty().SetSpecular(.4)
triangleEdgeActor.GetProperty().SetSpecularPower(10)
# Shrink the triangles we found earlier. Create the associated mapper
# and actor. Set the opacity of the shrunken triangles.
aShrinker = vtkShrinkPolyData()
aShrinker.SetShrinkFactor(1)
aShrinker.SetInputConnection(Marching.GetOutputPort())
aMapper = vtkPolyDataMapper()
aMapper.ScalarVisibilityOff()
aMapper.SetInputConnection(aShrinker.GetOutputPort())
Triangles = vtkActor()
Triangles.SetMapper(aMapper)
Triangles.GetProperty().SetDiffuseColor(
color.GetColor3d('banana'))
Triangles.GetProperty().SetOpacity(.6)
# Draw a cube the same size and at the same position as the one
# created previously. Extract the edges because we only want to see
# the outline of the cube. Pass the edges through a vtkTubeFilter so
# they are displayed as tubes rather than lines.
CubeModel = vtkCubeSource()
CubeModel.SetCenter(.5, .5, .5)
Edges = vtkExtractEdges()
Edges.SetInputConnection(CubeModel.GetOutputPort())
Tubes = vtkTubeFilter()
Tubes.SetInputConnection(Edges.GetOutputPort())
Tubes.SetRadius(.01)
Tubes.SetNumberOfSides(6)
Tubes.UseDefaultNormalOn()
Tubes.SetDefaultNormal(.577, .577, .577)
# Create the mapper and actor to display the cube edges.
TubeMapper = vtkPolyDataMapper()
TubeMapper.SetInputConnection(Tubes.GetOutputPort())
CubeEdges = vtkActor()
CubeEdges.SetMapper(TubeMapper)
CubeEdges.GetProperty().SetDiffuseColor(
color.GetColor3d('khaki'))
CubeEdges.GetProperty().SetSpecular(.4)
CubeEdges.GetProperty().SetSpecularPower(10)
# Create a sphere to use as a glyph source for vtkGlyph3D.
Sphere = vtkSphereSource()
Sphere.SetRadius(0.04)
Sphere.SetPhiResolution(20)
Sphere.SetThetaResolution(20)
# Remove the part of the cube with data values below 0.5.
ThresholdIn = vtkThresholdPoints()
ThresholdIn.SetInputData(Grid)
ThresholdIn.ThresholdByUpper(.5)
# Display spheres at the vertices remaining in the cube data set after
# it was passed through vtkThresholdPoints.
Vertices = vtkGlyph3D()
Vertices.SetInputConnection(ThresholdIn.GetOutputPort())
Vertices.SetSourceConnection(Sphere.GetOutputPort())
# Create a mapper and actor to display the glyphs.
SphereMapper = vtkPolyDataMapper()
SphereMapper.SetInputConnection(Vertices.GetOutputPort())
SphereMapper.ScalarVisibilityOff()
CubeVertices = vtkActor()
CubeVertices.SetMapper(SphereMapper)
CubeVertices.GetProperty().SetDiffuseColor(
color.GetColor3d('tomato'))
# Define the text for the label
caseLabel = vtkVectorText()
caseLabel.SetText('Case 1')
if label:
# Set up a transform to move the label to a new position.
aLabelTransform = vtkTransform()
aLabelTransform.Identity()
# Position the label according to the rotation of the figure.
if rotation == 0:
aLabelTransform.Translate(-0.2, 0, 1.25)
aLabelTransform.Scale(.05, .05, .05)
elif rotation == 1:
aLabelTransform.RotateY(90)
aLabelTransform.Translate(-1.25, 0, 1.25)
aLabelTransform.Scale(.05, .05, .05)
elif rotation == 2:
aLabelTransform.RotateY(180)
aLabelTransform.Translate(-1.25, 0, 0.2)
aLabelTransform.Scale(.05, .05, .05)
else:
aLabelTransform.RotateY(270)
aLabelTransform.Translate(-0.2, 0, 0.2)
aLabelTransform.Scale(.05, .05, .05)
# Move the label to a new position.
labelTransform = vtkTransformPolyDataFilter()
labelTransform.SetTransform(aLabelTransform)
labelTransform.SetInputConnection(caseLabel.GetOutputPort())
# Create a mapper and actor to display the text.
labelMapper = vtkPolyDataMapper()
labelMapper.SetInputConnection(labelTransform.GetOutputPort())
labelActor = vtkActor()
labelActor.SetMapper(labelMapper)
# Define the base that the cube sits on. Create its associated mapper
# and actor. Set the position of the actor.
baseModel = vtkCubeSource()
baseModel.SetXLength(1.5)
baseModel.SetYLength(.01)
baseModel.SetZLength(1.5)
baseMapper = vtkPolyDataMapper()
baseMapper.SetInputConnection(baseModel.GetOutputPort())
base = vtkActor()
base.SetMapper(baseMapper)
base.SetPosition(.5, -0.09, .5)
# Set the scalar values for this case of marching cubes.
# A negative case number will generate a complementary case
mcCase = mcCases[i]
if mcCase < 0:
cases[-mcCase](Scalars, caseLabel, 0, 1)
else:
cases[mcCase](Scalars, caseLabel, 1, 0)
# Force the grid to update.
Grid.Modified()
# Add the actors to the renderer
renderers[i].AddActor(triangleEdgeActor)
renderers[i].AddActor(base)
if label:
renderers[i].AddActor(labelActor)
renderers[i].AddActor(CubeEdges)
renderers[i].AddActor(CubeVertices)
renderers[i].AddActor(Triangles)
# Position the camera.
renderers[i].GetActiveCamera().Dolly(1.2)
# Rotate the camera an extra 30 degrees so the cube is not face on.
if rotation == 0:
renderers[i].GetActiveCamera().Azimuth(30)
elif rotation == 1:
renderers[i].GetActiveCamera().Azimuth(30 + 90)
elif rotation == 2:
renderers[i].GetActiveCamera().Azimuth(30 + 180)
else:
renderers[i].GetActiveCamera().Azimuth(30 + 270)
renderers[i].GetActiveCamera().Elevation(20)
renderers[i].ResetCamera()
renderers[i].ResetCameraClippingRange()
if i > 0:
renderers[i].SetActiveCamera(renderers[0].GetActiveCamera())
# Setup viewports for the renderers
rendererSize = 300
xGridDimensions = 4
if len(mcCases) < 4:
xGridDimensions = len(mcCases)
yGridDimensions = (len(mcCases) - 1) // 4 + 1
print('Grid dimensions, (x, y): ({:d}, {:d})'.format(xGridDimensions, yGridDimensions))
renWin.SetSize(
rendererSize * xGridDimensions, rendererSize * yGridDimensions)
renWin.SetWindowName('MarchingCases')
for row in range(0, yGridDimensions):
for col in range(0, xGridDimensions):
index = row * xGridDimensions + col
# (xmin, ymin, xmax, ymax)
viewport = [
float(col) / xGridDimensions,
float(yGridDimensions - (row + 1)) / yGridDimensions,
float(col + 1) / xGridDimensions,
float(yGridDimensions - row) / yGridDimensions]
renderers[index].SetViewport(viewport)
iren.Initialize()
renWin.Render()
iren.Start()
def case0(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, OUT)
scalars.InsertValue(1, OUT)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, OUT)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 0 - 00000000')
else:
caseLabel.SetText('Case 0c - 11111111')
def case1(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, OUT)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, OUT)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 1 - 00000001')
else:
caseLabel.SetText('Case 1c - 11111110')
def case2(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, IN)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, OUT)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 2 - 00000011')
else:
caseLabel.SetText('Case 2c - 11111100')
def case3(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, OUT)
scalars.InsertValue(2, IN)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, OUT)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 3 - 00000101')
else:
caseLabel.SetText('Case 3c - 11111010')
def case4(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, OUT)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, IN)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 4 - 01000001')
else:
caseLabel.SetText('Case 4c - 10111110')
def case5(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, OUT)
scalars.InsertValue(1, IN)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, IN)
scalars.InsertValue(5, IN)
scalars.InsertValue(6, OUT)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 5 - 00110010')
else:
caseLabel.SetText('Case 5c - 11001101')
def case6(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, OUT)
scalars.InsertValue(1, IN)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, IN)
scalars.InsertValue(4, IN)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, OUT)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 6 - 00011010')
else:
caseLabel.SetText('Case 6c - 11100101')
def case7(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, IN)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, IN)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 7 - 01000011')
else:
caseLabel.SetText('Case 7c - 10111100')
def case8(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, IN)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, IN)
scalars.InsertValue(5, IN)
scalars.InsertValue(6, OUT)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 8 - 00110011')
else:
caseLabel.SetText('Case 8c - 11001100')
def case9(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, OUT)
scalars.InsertValue(1, IN)
scalars.InsertValue(2, IN)
scalars.InsertValue(3, IN)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, IN)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 9 - 01001110')
else:
caseLabel.SetText('Case 9c - 10110001')
def case10(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, OUT)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, IN)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, IN)
scalars.InsertValue(6, IN)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 10 - 01101001')
else:
caseLabel.SetText('Case 10c - 10010110')
def case11(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, OUT)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, OUT)
scalars.InsertValue(4, IN)
scalars.InsertValue(5, IN)
scalars.InsertValue(6, IN)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 11 - 01110001')
else:
caseLabel.SetText('Case 11c - 10001110')
def case12(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, OUT)
scalars.InsertValue(1, IN)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, IN)
scalars.InsertValue(4, IN)
scalars.InsertValue(5, IN)
scalars.InsertValue(6, OUT)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 12 - 00111010')
else:
caseLabel.SetText('Case 12c - 11000101')
def case13(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, OUT)
scalars.InsertValue(1, IN)
scalars.InsertValue(2, OUT)
scalars.InsertValue(3, IN)
scalars.InsertValue(4, IN)
scalars.InsertValue(5, OUT)
scalars.InsertValue(6, IN)
scalars.InsertValue(7, OUT)
if IN == 1:
caseLabel.SetText('Case 13 - 01011010')
else:
caseLabel.SetText('Case 13c - 10100101')
def case14(scalars, caseLabel, IN, OUT):
scalars.InsertValue(0, IN)
scalars.InsertValue(1, OUT)
scalars.InsertValue(2, IN)
scalars.InsertValue(3, IN)
scalars.InsertValue(4, OUT)
scalars.InsertValue(5, IN)
scalars.InsertValue(6, IN)
scalars.InsertValue(7, IN)
if IN == 1:
caseLabel.SetText('Case 14 - 11101101')
else:
caseLabel.SetText('Case 14c - 00010010')
cases = [case0, case1, case2, case3, case4, case5, case6, case7, case8, case9, case10, case11, case12, case13, case14]
if __name__ == '__main__':
main()