1#!/usr/bin/env python
2import vtk
3from vtk.util.misc import vtkGetDataRoot
4VTK_DATA_ROOT = vtkGetDataRoot()
5
6# The resolution of the density function volume
7res = 100
8
9# Parameters for debugging
10NPts = 1000000
11math = vtk.vtkMath()
12math.RandomSeed(31415)
13
14# create pipeline
15#
16points = vtk.vtkBoundedPointSource()
17points.SetNumberOfPoints(NPts)
18points.ProduceRandomScalarsOn()
19points.ProduceCellOutputOff()
20points.Update()
21
22# Create a sphere implicit function
23sphere = vtk.vtkSphere()
24sphere.SetCenter(0.0,0.1,0.2)
25sphere.SetRadius(0.75)
26
27# Extract points within sphere
28extract = vtk.vtkFitImplicitFunction()
29extract.SetInputConnection(points.GetOutputPort())
30extract.SetImplicitFunction(sphere)
31extract.SetThreshold(0.005)
32extract.GenerateVerticesOn()
33
34# Clip out some of the points with a plane; requires vertices
35plane = vtk.vtkPlane()
36plane.SetOrigin(sphere.GetCenter())
37plane.SetNormal(1,1,1)
38
39clipper = vtk.vtkClipPolyData()
40clipper.SetInputConnection(extract.GetOutputPort())
41clipper.SetClipFunction(plane);
42
43# Generate density field from points
44# Use fixed radius
45dens0 = vtk.vtkPointDensityFilter()
46dens0.SetInputConnection(clipper.GetOutputPort())
47dens0.SetSampleDimensions(res,res,res)
48dens0.SetDensityEstimateToFixedRadius()
49dens0.SetRadius(0.05)
50dens0.SetDensityFormToVolumeNormalized()
51dens0.Update()
52vrange = dens0.GetOutput().GetScalarRange()
53
54map0 = vtk.vtkImageSliceMapper()
55map0.BorderOn()
56map0.SliceAtFocalPointOn()
57map0.SliceFacesCameraOn()
58map0.SetInputConnection(dens0.GetOutputPort())
59
60slice0 = vtk.vtkImageSlice()
61slice0.SetMapper(map0)
62slice0.GetProperty().SetColorWindow(vrange[1]-vrange[0])
63slice0.GetProperty().SetColorLevel(0.5*(vrange[0]+vrange[1]))
64
65# Now densify the point cloud and reprocess
66# Use relative radius
67print("Number of input points: {0}".format(clipper.GetOutput().GetNumberOfPoints()))
68denser = vtk.vtkDensifyPointCloudFilter()
69denser.SetInputConnection(clipper.GetOutputPort())
70denser.SetTargetDistance(0.025)
71denser.SetMaximumNumberOfIterations(5)
72denser.Update()
73print("Number of output points: {0}".format(denser.GetOutput().GetNumberOfPoints()))
74
75dens1 = vtk.vtkPointDensityFilter()
76dens1.SetInputConnection(denser.GetOutputPort())
77dens1.SetSampleDimensions(res,res,res)
78dens1.SetDensityEstimateToFixedRadius()
79dens1.SetRadius(0.05)
80dens1.SetDensityFormToVolumeNormalized()
81dens1.Update()
82vrange = dens1.GetOutput().GetScalarRange()
83
84map1 = vtk.vtkImageSliceMapper()
85map1.BorderOn()
86map1.SliceAtFocalPointOn()
87map1.SliceFacesCameraOn()
88map1.SetInputConnection(dens1.GetOutputPort())
89
90slice1 = vtk.vtkImageSlice()
91slice1.SetMapper(map1)
92slice1.GetProperty().SetColorWindow(vrange[1]-vrange[0])
93slice1.GetProperty().SetColorLevel(0.5*(vrange[0]+vrange[1]))
94
95# Create the RenderWindow, Renderer and both Actors
96#
97ren0 = vtk.vtkRenderer()
98ren0.SetViewport(0,0,0.5,1.0)
99ren1 = vtk.vtkRenderer()
100ren1.SetViewport(0.5,0,1,1.0)
101
102renWin = vtk.vtkRenderWindow()
103renWin.AddRenderer(ren0)
104renWin.AddRenderer(ren1)
105iren = vtk.vtkRenderWindowInteractor()
106iren.SetRenderWindow(renWin)
107
108# Add the actors to the renderer, set the background and size
109#
110ren0.AddActor(slice0)
111ren0.SetBackground(0,0,0)
112ren1.AddActor(slice1)
113ren1.SetBackground(0,0,0)
114
115renWin.SetSize(600,300)
116
117cam = ren0.GetActiveCamera()
118cam.ParallelProjectionOn()
119cam.SetFocalPoint(0,0,0)
120cam.SetPosition(0,0,1)
121ren0.ResetCamera()
122
123ren1.SetActiveCamera(cam)
124
125iren.Initialize()
126
127# render the image
128#
129renWin.Render()
130
131iren.Start()
132