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)
50#dens0.SetDensityEstimateToRelativeRadius()
51dens0.SetRelativeRadius(2.5)
52dens0.SetDensityFormToVolumeNormalized()
53
54# Time execution
55timer = vtk.vtkTimerLog()
56timer.StartTimer()
57dens0.Update()
58timer.StopTimer()
59time = timer.GetElapsedTime()
60print("Time to compute density field: {0}".format(time))
61vrange = dens0.GetOutput().GetScalarRange()
62print(dens0)
63
64map0 = vtk.vtkImageSliceMapper()
65map0.BorderOn()
66map0.SliceAtFocalPointOn()
67map0.SliceFacesCameraOn()
68map0.SetInputConnection(dens0.GetOutputPort())
69
70slice0 = vtk.vtkImageSlice()
71slice0.SetMapper(map0)
72slice0.GetProperty().SetColorWindow(vrange[1]-vrange[0])
73slice0.GetProperty().SetColorLevel(0.5*(vrange[0]+vrange[1]))
74
75# Generate density field from points
76# Use relative radius
77dens1 = vtk.vtkPointDensityFilter()
78dens1.SetInputConnection(clipper.GetOutputPort())
79dens1.SetSampleDimensions(res,res,res)
80#dens1.SetDensityEstimateToFixedRadius()
81dens1.SetRadius(0.05)
82dens1.SetDensityEstimateToRelativeRadius()
83dens1.SetRelativeRadius(2.5)
84dens1.SetDensityFormToNumberOfPoints()
85
86# Time execution
87timer = vtk.vtkTimerLog()
88timer.StartTimer()
89dens1.Update()
90timer.StopTimer()
91time = timer.GetElapsedTime()
92print("Time to compute density field: {0}".format(time))
93vrange = dens1.GetOutput().GetScalarRange()
94
95map1 = vtk.vtkImageSliceMapper()
96map1.BorderOn()
97map1.SliceAtFocalPointOn()
98map1.SliceFacesCameraOn()
99map1.SetInputConnection(dens1.GetOutputPort())
100
101slice1 = vtk.vtkImageSlice()
102slice1.SetMapper(map1)
103slice1.GetProperty().SetColorWindow(vrange[1]-vrange[0])
104slice1.GetProperty().SetColorLevel(0.5*(vrange[0]+vrange[1]))
105
106# Generate density field from points
107# Use fixed radius and weighted point density and volume normalized density
108# First need to generate some scalar attributes (weights)
109weights = vtk.vtkRandomAttributeGenerator()
110weights.SetInputConnection(clipper.GetOutputPort())
111weights.SetMinimumComponentValue(0.25)
112weights.SetMaximumComponentValue(1.75)
113weights.GenerateAllDataOff()
114weights.GeneratePointScalarsOn()
115
116dens2 = vtk.vtkPointDensityFilter()
117dens2.SetInputConnection(weights.GetOutputPort())
118dens2.SetSampleDimensions(res,res,res)
119dens2.SetDensityEstimateToFixedRadius()
120dens2.SetRadius(0.05)
121#dens2.SetDensityEstimateToRelativeRadius()
122dens2.SetRelativeRadius(2.5)
123dens2.SetDensityFormToVolumeNormalized()
124dens2.ScalarWeightingOn()
125
126# Time execution
127timer = vtk.vtkTimerLog()
128timer.StartTimer()
129dens2.Update()
130timer.StopTimer()
131time = timer.GetElapsedTime()
132print("Time to compute density field: {0}".format(time))
133vrange = dens2.GetOutput().GetScalarRange()
134
135map2 = vtk.vtkImageSliceMapper()
136map2.BorderOn()
137map2.SliceAtFocalPointOn()
138map2.SliceFacesCameraOn()
139map2.SetInputConnection(dens2.GetOutputPort())
140
141slice2 = vtk.vtkImageSlice()
142slice2.SetMapper(map2)
143slice2.GetProperty().SetColorWindow(vrange[1]-vrange[0])
144slice2.GetProperty().SetColorLevel(0.5*(vrange[0]+vrange[1]))
145
146# Generate density field from points
147# Use relative radius and weighted point density and npts density
148dens3 = vtk.vtkPointDensityFilter()
149dens3.SetInputConnection(weights.GetOutputPort())
150dens3.SetSampleDimensions(res,res,res)
151#dens3.SetDensityEstimateToFixedRadius()
152dens3.SetRadius(0.05)
153dens3.SetDensityEstimateToRelativeRadius()
154dens3.SetRelativeRadius(2.5)
155dens3.SetDensityFormToNumberOfPoints()
156dens3.ScalarWeightingOn()
157
158# Time execution
159timer = vtk.vtkTimerLog()
160timer.StartTimer()
161dens3.Update()
162timer.StopTimer()
163time = timer.GetElapsedTime()
164print("Time to compute density field: {0}".format(time))
165vrange = dens3.GetOutput().GetScalarRange()
166
167map3 = vtk.vtkImageSliceMapper()
168map3.BorderOn()
169map3.SliceAtFocalPointOn()
170map3.SliceFacesCameraOn()
171map3.SetInputConnection(dens3.GetOutputPort())
172
173slice3 = vtk.vtkImageSlice()
174slice3.SetMapper(map3)
175slice3.GetProperty().SetColorWindow(vrange[1]-vrange[0])
176slice3.GetProperty().SetColorLevel(0.5*(vrange[0]+vrange[1]))
177
178# Create the RenderWindow, Renderer and both Actors
179#
180ren0 = vtk.vtkRenderer()
181ren0.SetViewport(0,0,0.5,0.5)
182ren1 = vtk.vtkRenderer()
183ren1.SetViewport(0.5,0,1,0.5)
184ren2 = vtk.vtkRenderer()
185ren2.SetViewport(0,0.5,0.5,1)
186ren3 = vtk.vtkRenderer()
187ren3.SetViewport(0.5,0.5,1,1)
188
189renWin = vtk.vtkRenderWindow()
190renWin.AddRenderer(ren0)
191renWin.AddRenderer(ren1)
192renWin.AddRenderer(ren2)
193renWin.AddRenderer(ren3)
194iren = vtk.vtkRenderWindowInteractor()
195iren.SetRenderWindow(renWin)
196
197# Add the actors to the renderer, set the background and size
198#
199ren0.AddActor(slice0)
200ren0.SetBackground(0,0,0)
201ren1.AddActor(slice1)
202ren1.SetBackground(0,0,0)
203ren2.AddActor(slice2)
204ren2.SetBackground(0,0,0)
205ren3.AddActor(slice3)
206ren3.SetBackground(0,0,0)
207
208renWin.SetSize(300,300)
209
210cam = ren0.GetActiveCamera()
211cam.ParallelProjectionOn()
212cam.SetFocalPoint(0,0,0)
213cam.SetPosition(0,0,1)
214ren0.ResetCamera()
215
216ren1.SetActiveCamera(cam)
217ren2.SetActiveCamera(cam)
218ren3.SetActiveCamera(cam)
219
220iren.Initialize()
221
222# render the image
223#
224renWin.Render()
225
226iren.Start()
227