Visualizing streamlines using vtk and mevislab

30 Views Asked by At

I am working with a 4D CT scan of the lungs, and I am provided with 100 markers on the surface of the lungs. So basically, I have the coordinates of 100 points changing over 10 timesteps. I would like to plot streamlines for each of the points. I saw an example on GitHub, and I have changed it a bit. My code is mentioned below


import vtk
import numpy as np

# Is a nested list of length 10, and each element is a list of length 100
temporal_pts = []

NUMBER_OF_SEED_POINTS = 30

for i in range (10):
    filename = '/Markers/4DCT_Dicom_1_Points/' + str(i) + '0.pts'

    pts = []
    with open(filename, 'r') as file:
        for coordinate in file:
            coordinate = coordinate.strip()
            values = coordinate.split()
            coord = [float(values[0]), float(values[1]), float(values[2])]
            pts.append(coord)
        temporal_pts.append(pts)


points = vtk.vtkPoints()
source_points = vtk.vtkPoints()
iter = 0
for timesteps in temporal_pts:
    for pt in timesteps:
        points.InsertNextPoint(pt)
        if (iter < 100):
            source_points.InsertNextPoint(pt)
        iter += 1


# Vectors
vectors = vtk.vtkDoubleArray()
vectors.SetNumberOfComponents(3)
vectors.SetNumberOfTuples(points.GetNumberOfPoints())
for i in range(points.GetNumberOfPoints()):
    x, y, z = points.GetPoint(i)
    if (i+100 >= points.GetNumberOfPoints()):
        x1 = x
        y1 = y
        z1 = z
    else: 
        x1, y1, z1 = points.GetPoint(i+100)
    vectors.SetTuple3(i, (x1-x), (y1-y), (z1-z) )  



# Create a structured grid
grid = vtk.vtkStructuredGrid()
grid.SetDimensions(21, 21, 21)
grid.SetPoints(points)
grid.GetPointData().SetVectors(vectors)


# Create a polydata object
source_polydata = vtk.vtkPolyData()
source_polydata.SetPoints(source_points)

# Create a source using the polydata
seeds = vtk.vtkPointSource()
seeds.SetCenter(0, 0, 0)
seeds.SetNumberOfPoints(source_points.GetNumberOfPoints())  # Set the number of points
seeds.SetRadius(1.0)  # Set the radius to 1 or any other value
seeds.Update()


streamer = vtk.vtkStreamTracer()
streamer.SetInputData(grid)
streamer.SetSourceConnection(seeds.GetOutputPort())
streamer.SetMaximumPropagation(500)  # Increase this to let the streamlines go further.
streamer.SetIntegrationStepUnit(vtk.vtkStreamTracer.LENGTH_UNIT)
streamer.SetInitialIntegrationStep(0.1)  # Increased step size
streamer.SetIntegrationDirectionToBoth()
streamer.SetComputeVorticity(False)

# Visualize the streamlines.
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(streamer.GetOutputPort())
actor = vtk.vtkActor()
actor.SetMapper(mapper)

renderer = vtk.vtkRenderer()
renderer.AddActor(actor)
renderer.SetBackground(0, 0, 0)

# Create a glyph source (a sphere)
glyph_source = vtk.vtkSphereSource()
glyph_source.SetRadius(0.5)

# Create a polydata object
source_polydata = vtk.vtkPolyData()
source_polydata.SetPoints(source_points)

render_window = vtk.vtkRenderWindow()
render_window.AddRenderer(renderer)

interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(render_window)
interactor.Initialize()
interactor.Start()

The vector for each point will be the coordinates of that point in the next timestep minus the coordinates at that timestep (The vector for the last timestep would be zero) .But, the streamlines come out looking weird, and I have looked at the documentation but couldnt find a solution. I think I am going a bit wrong in code after creating the vectors. Any kind of help would be appreciated. Thanks

0

There are 0 best solutions below