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()