How to volume render 3D numpy array using VTK (MarchingCubes) in Python?

Question

I have a 3D numpy array and I am trying to volume render it using VTK. However, I get a completely different volume rendering when I visualise it. I suspect it has something to do with my conversion of numpy array to the VTK image format but I can't seem to figure out where I am going wrong. I have uploaded the numpy array here.

Can someone help me figure out where I'm going wrong?

This is my code:

#!/usr/bin/env python

import os
import numpy as np

ArrayDicom = np.load('test3.npy')
data_matrix = ArrayDicom
w, d, h = ArrayDicom.shape

colors = vtkNamedColors()
iso_value = 200

reader = vtkImageImport()
data_string = data_matrix.tobytes()
reader.CopyImportVoidPointer(data_string, len(data_string))
reader.SetDataScalarTypeToUnsignedChar()
reader.SetNumberOfScalarComponents(1)
reader.SetDataExtent(0, w-1, 0, d-1, 0, h-1)
reader.SetWholeExtent(0, w-1, 0, d-1, 0, h-1)
reader.Update()

volume = vtkImageData()
volume.DeepCopy(reader.GetOutput())

surface = vtkMarchingCubes()
surface.SetInputData(volume)
surface.ComputeNormalsOn()
surface.SetValue(0, iso_value)

renderer = vtkRenderer()
renderer.SetBackground(colors.GetColor3d('DarkSlateGray'))

render_window = vtkRenderWindow()
render_window.AddRenderer(renderer)
render_window.SetWindowName('MarchingCubes')

interactor = vtkRenderWindowInteractor()
interactor.SetRenderWindow(render_window)

mapper = vtkPolyDataMapper()
mapper.SetInputConnection(surface.GetOutputPort())
mapper.ScalarVisibilityOff()

actor = vtkActor()
actor.SetMapper(mapper)
actor.GetProperty().SetColor(colors.GetColor3d('MistyRose'))

renderer.AddActor(actor)

render_window.Render()
interactor.Start()

This is my volume rendering:

This is my expected volume rendering:

Dave Chen · Accepted Answer · 2021-11-12 17:32:08Z

1

Numpy uses a different array ordering than VTK. You should be able to re-order w, h and d to get the right thing.

This is how you want it:

h, d, w = ArrayDicom.shape

OK, here's a conversion script that I used to convert to a VTK file:

import numpy as np
import SimpleITK as sitk

x = np.load("test3.npy")
y = sitk.GetImageFromArray(x)
sitk.WriteImage(y, "test3.vtk")

It's not as nice as correctly getting the VTK image import to work, but, well, I'm a SimpleITK guy, and I know that converting numpy works in SimpleITK.

edited Nov 12, 2021 at 17:32

answered Nov 12, 2021 at 16:32

Dave Chen

2,1251 gold badge14 silver badges21 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

user42 Over a year ago

Hi Dave, thank you for your answer. I tried changing h, d and w but it doesn't seem to work, I still get a volume with the same square shape as in the image. Could it be due to the SetDataExtent and SetWholeExtent or the spacing?

user42 Over a year ago

Hi Dave, I will try that as well, thank you once again.

Collectives™ on Stack Overflow

How to volume render 3D numpy array using VTK (MarchingCubes) in Python?

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related