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