Calculate volumes counting voxels, is a correct approach?

Hello,

I have to calculate the volume of some lesion masks and I was wondering if counting the voxels != 0 and multiplying for the dimension of a single voxel in mm³ can be a right approach.

I mean:

  • Load the Nifti with nibabel
  • Iterating x,y,z and counting the voxel != 0
  • Multiplying for pixmap[1] * pixmap[2] + pixmap[3] from Nifti header
1 Like

@aorco you could do something like:

import nibabel as nb
import numpy as np
img = nb.load("path/to/img.nii")
num_vox = np.sum(img.get_data() > 0)
2 Likes

In the case you have FSL installed, you can use fslstats to compute that for you.

The usage syntax would be:

$ fslstats input.nii.gz -V

Note the capital V. This will output two numbers, the first is the number of voxels greater than 0, the second is their volume in cubic millimiters. Ex:

$fslstats $FSLDIR/data/standard/LowerCingulum_1mm.nii.gz -V
10964 10964.166992

You can check the pixel dimensions of the file to see that the volume calculation is correct.

2 Likes

Is there already a python implementation / wrapper for this?
Would be happy to write one otherwise. In which package would this fit best?

If you’re making a general equivalent to fslstats or fslmaths, then nibabel would be a reasonable place for these. We have a number of CLI tools following the pattern nib-*.

for now I was thinking of simply reimplementing fslstats -V
will have a look where it could fit in

added a draft PR to nibabel:

1 Like

Number of voxels are only equal to the volume when the voxels are 1 mm^3. For smaller/larger and/or anisotropic voxels, something like this would work when using nibabel:

"""Compute volume of non-zero voxels in a nifti image."""

import numpy as np
import nibabel as nb

INPUT = "/path/to/image.nii.gz"

# Load data
nii = nb.load(INPUT)
img = nii.get_fdata()

# Get voxel dimensions
voxel_dims = (nii.header["pixdim"])[1:4]
print("Voxel dimensions:")
print("  x = {} mm".format(voxel_dims[0]))
print("  y = {} mm".format(voxel_dims[1]))
print("  z = {} mm".format(voxel_dims[2]))

# Compute volume
nonzero_voxel_count = np.count_nonzero(img)
voxel_volume = np.prod(voxel_dims)
nonzero_voxel_volume = nonzero_voxel_count * voxel_volume

print("Number of non-zero voxels = {}".format(nonzero_voxel_count))
print("Volume of non-zero voxels = {} mm^3".format(nonzero_voxel_volume))

Thank you very much! it works :slightly_smiling_face:

Just to make it more robust, I think it would be safer to change

nonzero_voxel_count = np.sum(img > 0)

with
nonzero_voxel_count = np.count_nonzero(img)

because the segmentation label can also be a value different than 1, and in that case sum can produce misleading results.

1 Like

Thanks @Tommaso_Di_Noto , you are right. Updated my post :slight_smile:

Edit: Actually I think you were wrong because img > 0 returns booleans so segmentation labels higher than 1 become True which is 1 anyway when summing. However, I would be missing the negative values so np.count_nonzero(img) or np.sum(img != 0) would still improve upon np.sum(img > 0) :slight_smile:

Ah you are right! Even better like this @ofgulban :slight_smile:

1 Like

Maybe a little late for this, but 3dBrickStat in AFNI can provide either number of voxels or volume.

https://afni.nimh.nih.gov/afni/community/board/read.php?1,77807,77820#msg-77820

MANAGED BY INCF