Conversion of NIfTI mask file into VTK polydata file

sunita · July 10, 2023, 10:31am

Dear VTK users,

We are trying to convert 3D numpy array into VTK polydata file to be used in three.js vtk loader which only supports polydata file.

When we generated vtk file using ITK-SNAP we could load it to vtk loader. However, when we used following code we generated .vtk file but it is not loading to vtk loader and we are getting unsupported dataset type structured points. What could be the reason? Where we are going wrong?

def convert_nifti_to_vtk(nifti_file):
reader = vtk.vtkNIFTIImageReader()
reader.SetFileName(nifti_file)
reader.Update()

polydata = vtk.vtkPolyData()
vtk.vtkMarchingCubes().SetInputData(reader.GetOutput())
vtk.vtkMarchingCubes().Update()
polydata.ShallowCopy(vtk.vtkMarchingCubes().GetOutput())

return polydata

if name == “main”:
nifti_file = “my_nifti_file.nii”
polydata = convert_nifti_to_vtk(nifti_file)

writer = vtk.vtkPolyDataWriter()
writer.SetFileName(“my_vtk_file.vtk”)
writer.SetInputData(polydata)
writer.Write()

dgobbi · July 10, 2023, 7:37pm

Looks like a typo. You are calling the vtkMarchingCubes constructor several times

vtk.vtkMarchingCubes().SetInputData(reader.GetOutput())
vtk.vtkMarchingCubes().Update()
polydata.ShallowCopy(vtk.vtkMarchingCubes().GetOutput())

The constructor should be called just once, to create the object:

mc = vtk.vtkMarchingCubes()
mc.SetInputData(reader.GetOutput())
mc.SetValue(0, 0.5) # isocontour value
mc.Update()

Also note that NIFTI images have transform matrices that should be applied to the data

transform = vtk.vtkTransform()
if reader.GetQFormMatrix():
    transform.Concatenate(reader.GetQFormMatrix())
elif reader.GetSFormMatrix():
    transform.Concatenate(reader.GetSFormMatrix())

tpd = vtk.vtkTransformPolyData()
tpd.SetInputData(mc.GetOutput())
tpd.SetTransform(transform)
tpd.Update()

polydata.ShallowCopy(tpd.GetOutput())

sunita · July 12, 2023, 12:54pm

Hi David,

Thanks for your help. We updated the code as per suggestions. We could generate vtk file now and load it into three.js vtkloader. However, the edges are not smooth. We are getting imperfect image because of the rugged edges. Then, we used smoothing functions that come with vtk and generated the .vtk file. But image is not loading to three.js vtkloader.

Is the code correct or we are missing something?

Reader

reader = vtk.vtkNIFTIImageReader()
reader.SetFileName(‘./nf_ed.nii.gz’)
reader.Update()

Marching Cubes

mc=vtk.vtkMarchingCubes()
mc.SetInputData(reader.GetOutput())
mc.ComputeNormalsOn()
mc.ComputeGradientsOn()
mc.SetValue(0,0.5)
mc.Update()

Transform

transform = vtk.vtkTransform()
if reader.GetQFormMatrix():
transform.Concatenate(reader.GetQFormMatrix())
elif reader.GetSFormMatrix():
transform.Concatenate(reader.GetSFormMatrix())

tpd = vtk.vtkTransformPolyDataFilter()
tpd.SetInputData(mc.GetOutput())
tpd.SetTransform(transform)
tpd.Update()

Windowed SincPolyData Filter

smt = vtk.vtkWindowedSincPolyDataFilter()
smt.SetInputConnection(tpd.GetOutputPort())
smt.SetNumberOfIterations(20)
smt.SetPassBand(0.1)
smt.NormalizeCoordinatesOn()
smt.FeatureEdgeSmoothingOn()
smt.SetFeatureAngle(120.0)
smt.SetEdgeAngle(90.0)
smt.BoundarySmoothingOn()
smt.NonManifoldSmoothingOn()

PolyData Mapper

mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(smt.GetOutputPort())
mapper.ScalarVisibilityOff()

PolyData Writer

writer = vtk.vtkPolyDataWriter()
writer.SetInputData(mapper.GetInput())
writer.SetFileName(‘./new_mask.vtk’)
writer.Update()

Kindly please guide.

Thank you David.

Best regards,
Sunita

dgobbi · July 12, 2023, 1:55pm

Hi Sunita,

The best option is to smooth the mask image before applying vtkMarchingCubes. This generally gives better results than smoothing the polydata afterwards.

So I would suggest converting the image to floating-point, smoothing it, and then contouring it with an isovalue that is exactly halfway between the min and max values in the mask. Some masks have a range of [0,1], while others have a range of [0,255].

cast = vtk.vtkImageCast()
cast.SetInputData(reader.GetOutput())
cast.SetOutputScalarTypeToFloat()
cast.Update()

smooth = vtk.vtkImageGaussianSmooth()
smooth.SetInputData(cast.GetOutput())
smooth.SetStandardDeviations(1.2, 1.2, 1.2) # see note 1
smooth.Update()

mc=vtk.vtkMarchingCubes()
mc.SetInputData(smooth.GetOutput())
mc.SetValue(0,0.5) # see note 2
mc.Update()

Standard deviation is in voxel units, e.g. this sets the standard deviation to 1.2 times the voxel size
Set the isovalue according to the range of the mask image: 0.5 for [0,1], or 127.5 for [0,255].

Regarding the problem with the saved data: it’s due to this code:

writer.SetInputData(mapper.GetInput())

You haven’t called Update() on your vtkWindowedSincPolyDataFilter, so the writer is writing the data before the data is ready. There is information on this discussion list and elsewhere on the differences between SetInputData() vs. SetInputConnection() and how they effect VTK pipeline execution.

Edit: fixed Note 1 regarding units of StandardDeviation

sunita · July 13, 2023, 8:32am

Hi David,

Thanks a lot.

We could generate a smooth image of the mask following your suggestion. It’s working great. The code we are using now is given below.

Reader

reader = vtk.vtkNIFTIImageReader()
reader.SetFileName(msk_path)
reader.Update()

Image Cast

cast = vtk.vtkImageCast()
cast.SetInputData(reader.GetOutput())
cast.SetOutputScalarTypeToFloat()
cast.Update()

Smooth

smooth = vtk.vtkImageGaussianSmooth()
smooth.SetInputData(cast.GetOutput())
smooth.SetStandardDeviations(1.2, 1.2, 1.2) # see note 1
smooth.Update()

Marching Cubes

mc=vtk.vtkMarchingCubes()
mc.SetInputData(smooth.GetOutput())
mc.ComputeNormalsOn()
mc.SetValue(0,0.5)
mc.Update()

Transform

transform = vtk.vtkTransform()
if reader.GetQFormMatrix():
    transform.Concatenate(reader.GetQFormMatrix())
elif reader.GetSFormMatrix():
    transform.Concatenate(reader.GetSFormMatrix())

tpd = vtk.vtkTransformPolyDataFilter()
tpd.SetInputData(mc.GetOutput())
tpd.SetTransform(transform)
tpd.Update()

PolyData Writer

writer = vtk.vtkPolyDataWriter()
writer.SetInputData(tpd.GetOutput())
writer.SetFileName(msk_name)
writer.Update()

But still there is a problem. We are unable to load the smooth image into three.js. It looks same as earlier which has rough surface as shown in the attached file.

The code we are using to load is as follows.

import * as THREE from 'three';
import { TrackballControls } from 'three/addons/controls/TrackballControls.js';
import { VTKLoader } from 'three/ addons/loaders/VTKLoader.js';

const loader_1 = new VTKLoader();
loader_1.load( './file_ed.vtk', function ( geometry ) {
	geometry.center();
	geometry.computeVertexNormals();
	const material = new THREE.MeshLambertMaterial( { color: 0x00ff00 } );
	const mesh = new THREE.Mesh( geometry, material );
	scene.add( mesh );
});

We will highly appreciate your help David.

Thanks again for your prompt reply.

Best regards,
Sunita

dgobbi · July 13, 2023, 12:18pm

Since the rendered data is colored, my guess is that the .vtk file has scalar values associated with the points in the dataset.

Try adding the following line of code for marching cubes, so that it doesn’t add scalars to the data:

mc.ComputeScalarsOff()

Edit: It looks like you are setting the color to green ( { color: 0x00ff00 } ) so maybe the scalars aren’t the problem. But try turning them off anyway, just in case.

sunita · July 13, 2023, 1:02pm

Hi David,

We added the line. But still it is the same.
We are unable to load the smooth vtk file to three.js vtkloader.

We removed the green color and now the default color is white but the image looks same.

Please suggest where it is going wrong or any other aspect we need to look into.

Thank you so much.

Best regards,
Sunita

dgobbi · July 13, 2023, 1:25pm

You can also try turning off normals:

mc.ComputeNormalsOff()

I don’t have any other suggestions. It looks like a rendering problem in three.js, but I have never used three.js.

Sebastien_Jourdain · July 13, 2023, 5:27pm

Did you try to load it here to make sure your vtp file is as expected?

sunita · July 13, 2023, 5:39pm

You are right David. The problem is with three.js. A vtk file obtained from Simple ITK also shows rough surface when loaded using three.js. We need to figure out the problem.

Just a few more questions. How to make the image semitransparent using vtk? Can we use vtk libraries to make mesh plot or dot plot?

Your guidance helped a lot to come to this point. Before we asked the question in the vtk forum, we searched a lot but didn’t find much help.

Thank you so very much and wishing you all the best for your future endeavor.

Best regards,
Sunita

sunita · July 14, 2023, 7:53am

Hi David,

Finally, we could dig into the problem. The problem is with the vtk version. vtk file generated with older version that is 6.3.0 can be loaded into three.js without any problem. However, when it is created using vtk 9.1.0 version or higher it is showing discontinuous and rough surface.

We created two versions of vtk file using different versions of ITK-SNAP. ITK-SNAP 3.8.0 which uses vtk 6.3.0 generates the the vtk image which can be loaded into three.js vtk loader and it appears smooth and continuous. While the vtk image generated using ITK-SNAP 4.0.1 which uses vtk 9.1.0 shows discontinuous and rough surface.

Could you please help us in getting lower version of vtk that is vtk 6.3.0 ?

We are stuck in this problem quite sometime now. Please help.

We are very close to solve this problem. However, we can not do unless we will get older version of it.

Thanks a lot.

Best regards,
Sunita

sunita · July 14, 2023, 9:14am

Hi Sebastien,

The vtk files generated using 6.3.0 or 9.1.0 can be seen in paraview without any problem. It is the three.js which is not supporting vtk file generated with 9.1.0 version.

Thanks for the link. vtk converted to vtp can seen in that link which shows smooth rendering.

May be we have to implement vtk.js.

Thank you so much.

Sunita

dgobbi · July 14, 2023, 11:38am

Hi Sunita,

An old version of VTK is not needed, because the writer can create the old VTK file format. For example,

writer.SetFileVersion(writer.VTK_LEGACY_READER_VERSION_5_1)
writer.Write()

or

writer.SetFileVersion(writer.VTK_LEGACY_READER_VERSION_4_2)
writer.Write()

One of these should work.

Regarding your other questions, yes, it is possible to make an image semitransparent in VTK either by using the alpha channel (RGBA rendering), or by calling SetOpacity(alpha) on the image property, using an alpha between 0 and 1.

VTK provides multiple ways of doing mesh and dot plots, depending on the requirements. The vtkAxisActor and vtkCubeAxesActor are also useful for showing gridlines for 3D plots.

sunita · July 15, 2023, 6:51am

Hi Sebastien,

Thanks for your suggestion.

We are now trying to implement vtk.js using CDN. Actor and mapper are imported successfully. Which reader is used to read vtp file and how to import using CDN?

Your help will be highly appreciated.

Best regards,
Sunita