Different Volume Calculation Results

Development
1

I have a binary mask, I loaded it in 3D slicer as a segmentation file, and in Data module, I used “export visible segments to model”. Then in Models module, I can see the model volume is 55349.06 mm^3 (none smooth).

image
image1348×392 105 KB

Then, I used my own vtk codes to calculate its volume:

mask = vtk.vtkNIFTIImageReader()
mask.SetFileName('origCT_nodule_0.nii.gz')
mask.Update()
nodule = vtk.vtkDiscreteMarchingCubes()
nodule.SetInputData(mask.GetOutput())
nodule.Update()

massProperties = vtk.vtkMassProperties()
massProperties.SetInputData(nodule.GetOutput())
massProperties.Update()
volume = massProperties.GetVolume()
print(volume)

However, the volume is 47789.92 mm^3. I visualized it to check what happened:

image
image1135×701 65.7 KB

Oh, it seems this is an unclosed surface representation. I know 3D Slicer will close it first and then calculate its volume (as in the first screenshot). Then, I used vtkFillHolesFilter() to close the surface first:

fill_holes_filter = vtk.vtkFillHolesFilter()
fill_holes_filter.SetInputData(nodule.GetOutput())
fill_holes_filter.SetHoleSize(100.0)
fill_holes_filter.Update()

image
image1135×701 50.4 KB

This time, the volume is 40265.35 mm^3. What happened???

2

This class needs also a triangulated surface, try applying vtkTriangleFilter().

3

Estimating volume from a mesh is much more complicated and many things can influence the result. Slicer uses more robust and accurate surface reconstruction method than discrete marching cubes. If you use a simpler method then you will get a different, most likely less accurate results.

If you have a labelmap image then you can compute the volume very easily and robustly, using Segment Statistics module. There is only one single computation method (multiply the volume of one voxel by the number of voxels in the segment), so there are no options, no questions, no doubts.

4

I tried the following codes but there’s no luck

t = vtk.vtkTriangleFilter()
t.SetInputData(nodule.GetOutput())
t.PassVertsOn()
t.PassLinesOn()
t.Update()

5

Solved it by using the following codes, thanks for the idea:

image_data = mask_image2.GetOutput()
spacing = image_data.GetSpacing()

# Extract the voxel data and convert to a NumPy array
dims = image_data.GetDimensions()
scalars = image_data.GetPointData().GetScalars()
from vtkmodules.util import numpy_support
mask_image_array = numpy_support.vtk_to_numpy(scalars)
mask_image_array = mask_image_array.reshape(dims, order='F')  # Reshape to 3D array with Fortran order

# Convert the segmentation mask to a binary array (1 for segmented regions, 0 otherwise)
binary_array = np.where(mask_image_array > 0, 1, 0)

# Count the number of non-zero voxels
voxel_count = binary_array.sum()

# Calculate the volume in physical units (e.g., mm^3)
voxel_volume = spacing[0] * spacing[1] * spacing[2]
total_volume = voxel_count * voxel_volume

print(total_volume)