KochSnowflake
VTKExamples/Python/DataManipulation/KochSnowflake
Description
==Description== This script draws a Koch snowflake using the VTK. The general idea is to exercise some of the components of a vtkPolyData to produce something interesting rather than a boring old cube. Not that I have anything against cubes. There is also a C++ version of this example: http://www.vtk.org/Wiki/VTK/Examples/Cxx/Visualization/KochSnowflake
Code
KochSnowflake.py
#!/usr/bin/env python
#------------------------------------------------------------------------------#
# Imports #
#------------------------------------------------------------------------------#
from math import pi, cos, sin, sqrt
import vtk
LEVEL = 6
#------------------------------------------------------------------------------#
# Koch Snowflake as vtkPolyLine #
#------------------------------------------------------------------------------#
def as_polyline(points, level):
# Use the points from the previous iteration to create the points of the next
# level. There is an assumption on my part that the curve is traversed in a
# counterclockwise fashion. If the initial triangle above is written to
# describe clockwise motion, the points will face inward instead of outward.
for i in range(level):
temp = vtk.vtkPoints()
# The first point of the previous vtkPoints is the first point of the next vtkPoints.
temp.InsertNextPoint(*points.GetPoint(0))
# Iterate over "edges" in the vtkPoints
for i in range(1, points.GetNumberOfPoints()):
x0, y0, z0 = points.GetPoint(i - 1)
x1, y1, z1 = points.GetPoint(i)
t = sqrt((x1 - x0)**2 + (y1 - y0)**2)
nx = (x1 - x0)/t # x-component of edge unit tangent
ny = (y1 - y0)/t # y-component of edge unit tangent
# the points describing the Koch snowflake edge
temp.InsertNextPoint(x0 + nx*t/3, y0 + ny*t/3, 0.)
temp.InsertNextPoint(x0 + nx*t/2 + ny*t*sqrt(3)/6, y0 + ny*t/2 - nx*t*sqrt(3)/6, 0.)
temp.InsertNextPoint(x0 + nx*2*t/3, y0 + ny*2*t/3, 0.)
temp.InsertNextPoint(x0 + nx*t, y0 + ny*t, 0.)
points = temp
# draw the outline
lines = vtk.vtkCellArray()
pl = vtk.vtkPolyLine()
pl.GetPointIds().SetNumberOfIds(points.GetNumberOfPoints())
for i in range(points.GetNumberOfPoints()):
pl.GetPointIds().SetId(i, i)
lines.InsertNextCell(pl)
# complete the polydata
polydata = vtk.vtkPolyData()
polydata.SetLines(lines)
polydata.SetPoints(points)
return polydata
#------------------------------------------------------------------------------#
# Koch Snowflake as collection of vtkTriangles #
#------------------------------------------------------------------------------#
def as_triangles(indices, cellarray, level, data):
if len(indices) >= 3:
stride = len(indices)//4
indices.append(indices[-1] + 1)
triangle = vtk.vtkTriangle()
triangle.GetPointIds().SetId(0, indices[ stride])
triangle.GetPointIds().SetId(1, indices[2*stride])
triangle.GetPointIds().SetId(2, indices[3*stride])
cellarray.InsertNextCell(triangle)
data.InsertNextValue(level)
as_triangles(indices[ 0: stride], cellarray, level + 1, data)
as_triangles(indices[ stride: 2*stride], cellarray, level + 1, data)
as_triangles(indices[2*stride: 3*stride], cellarray, level + 1, data)
as_triangles(indices[3*stride: -1], cellarray, level + 1, data)
#------------------------------------------------------------------------------#
# Main Method #
#------------------------------------------------------------------------------#
if __name__ == "__main__":
# Initially, set up the points to be an equilateral triangle. Note that the
# first point is the same as the last point to make this a closed curve when
# I create the vtkPolyLine.
points = vtk.vtkPoints()
for i in range(4):
points.InsertNextPoint(cos(2.*pi*i/3), sin(2*pi*i/3.), 0.)
outline_pd = as_polyline(points, LEVEL)
# You have already gone through the trouble of putting the points in the
# right places - so "all" you need todo now is to create polygons from the
# points that are in the vtkPoints.
# The points that are passed in, have an overlap of the beginning and the
# end. For this next trick, I will need a list of the indices in the
# vtkPoints. They're consecutive, so thats pretty straightforward.
indices = [i for i in range(outline_pd.GetPoints().GetNumberOfPoints() + 1)]
triangles = vtk.vtkCellArray()
# Set this up for each of the initial sides, then call the recursive function.
stride = (len(indices) - 1)//3
# The cell data will allow us to color the triangles based on the level of
# the iteration of the Koch snowflake.
data = vtk.vtkIntArray()
data.SetNumberOfComponents(0)
data.SetName("Iteration Level")
# This is the starting triangle.
t = vtk.vtkTriangle()
t.GetPointIds().SetId(0, 0)
t.GetPointIds().SetId(1, stride)
t.GetPointIds().SetId(2, 2*stride)
triangles.InsertNextCell(t)
data.InsertNextValue(0)
as_triangles(indices[ 0: stride + 1], triangles, 1, data)
as_triangles(indices[ stride: 2*stride + 1], triangles, 1, data)
as_triangles(indices[2*stride: -1], triangles, 1, data)
triangle_pd = vtk.vtkPolyData()
triangle_pd.SetPoints(outline_pd.GetPoints())
triangle_pd.SetPolys(triangles)
triangle_pd.GetCellData().SetScalars(data)
#-----------------#
# rendering stuff #
#-----------------#
outline_mapper = vtk.vtkPolyDataMapper()
outline_mapper.SetInputData(outline_pd)
lut = vtk.vtkLookupTable()
lut.SetNumberOfTableValues(256)
lut.SetHueRange(0.6, 0.6)
lut.SetSaturationRange(0., 1.)
lut.Build()
triangle_mapper = vtk.vtkPolyDataMapper()
triangle_mapper.SetInputData(triangle_pd)
triangle_mapper.SetScalarRange(0.0, LEVEL)
triangle_mapper.SetLookupTable(lut)
outline_actor = vtk.vtkActor()
outline_actor.SetMapper(outline_mapper)
triangle_actor = vtk.vtkActor()
triangle_actor.SetMapper(triangle_mapper)
outline_ren = vtk.vtkRenderer()
outline_ren.AddActor(outline_actor)
outline_ren.SetViewport(0.0, 0.0, 0.5, 1.0)
triangle_ren = vtk.vtkRenderer()
triangle_ren.AddActor(triangle_actor)
triangle_ren.SetViewport(0.5, 0.0, 1.0, 1.0)
triangle_ren.SetActiveCamera(outline_ren.GetActiveCamera())
renw = vtk.vtkRenderWindow()
renw.AddRenderer(outline_ren)
renw.AddRenderer(triangle_ren)
renw.SetSize(800, 400)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renw)
outline_ren.ResetCamera()
renw.Render()
iren.Start()