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