Good algorithm for computing volume or surface area in python

Question

I am trying to compute the volume (or surface area) of a 3D numpy array. The voxels are anisotropic in a lot of cases, and I have the pixel to cm conversion factor in each direction.

Does anyone know a good place to find a toolkit or package to do the above??

Right now, I have some in-house code, but I am looking to upgrade to something more industrial strength in terms of accuracy.

Edit1: Here is some (poor) sample data. This is much smaller than the typical sphere. I will add better data when I can generate it! It is in (self.)tumorBrain.tumor.

You need to define what you mean a bit more clearly. If you data is actually a 3D array, then the volume the entire grid occupies is nx * dx * ny * dy * nz * dz, but I'm pretty sure you didn't mean that... Do you want the volume inside an isosurface? — Joe Kington
– Joe Kington, Commented Jan 18, 2012 at 0:56
I THINK that you are correct. It is a X x Y x Z binary array and I want the volume of everything contained in the perimeter of the binary part of it. It is typically (but not always) sphere shaped. — tylerthemiler
– tylerthemiler, Commented Jan 18, 2012 at 1:01
this sounds like fun, can you post a link to some example data? just save the numpy array with pickle.dump — wim
– wim, Commented Jan 18, 2012 at 1:14
@tylerthemiller - Just FYI: If you pickle numpy arrays, be sure to specify specify something other than the default ascii protocol. Otherwise, you'll wind up with a huge file (like your example). E.g. pickle.dump(data, f, pickle.HIGHEST_PROTOCOL) Alternately, you can just use numpy.save. Also, we can't load your data because you pickled some sort of dataStructures object, and without your code, we can't unpickle it. Try pickling just the numpy array and a tuple of (dx, dy, dz). — Joe Kington
– Joe Kington, Commented Jan 18, 2012 at 2:10
Since you are dealign with medical images, you should take a look at invesalius - A Free software for 3d medical imaging written in Python (sponsored by the brazilian government) - it will likely hae the features you need, or allow you to code then in. softwarepublico.gov.br/ver-comunidade?community_id=626732 — jsbueno
– jsbueno, Commented Jan 18, 2012 at 2:55

Joe Kington · Accepted Answer · 2012-01-18 14:41:38Z

5

One option is to use VTK. (I'm going to use the tvtk python bindings for it here...)

At least in some circumstances, getting the area within the isosurface will be a bit more accurate.

Also, as far as surface area goes, tvtk.MassProperties calculates surface area as well. It's mass.surface_area (with the mass object in the code below).

import numpy as np
from tvtk.api import tvtk

def main():
    # Generate some data with anisotropic cells...
    # x,y,and z will range from -2 to 2, but with a 
    # different (20, 15, and 5 for x, y, and z) number of steps
    x,y,z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j]
    r = np.sqrt(x**2 + y**2 + z**2)

    dx, dy, dz = [np.diff(it, axis=a)[0,0,0] for it, a in zip((x,y,z),(0,1,2))]

    # Your actual data is a binary (logical) array
    max_radius = 1.5
    data = (r <= max_radius).astype(np.int8)

    ideal_volume = 4.0 / 3 * max_radius**3 * np.pi
    coarse_volume = data.sum() * dx * dy * dz
    est_volume = vtk_volume(data, (dx, dy, dz), (x.min(), y.min(), z.min()))

    coarse_error = 100 * (coarse_volume - ideal_volume) / ideal_volume
    vtk_error = 100 * (est_volume - ideal_volume) / ideal_volume

    print 'Ideal volume', ideal_volume
    print 'Coarse approximation', coarse_volume, 'Error', coarse_error, '%'
    print 'VTK approximation', est_volume, 'Error', vtk_error, '%'

def vtk_volume(data, spacing=(1,1,1), origin=(0,0,0)):
    data[data == 0] = -1
    grid = tvtk.ImageData(spacing=spacing, origin=origin)
    grid.point_data.scalars = data.T.ravel() # It wants fortran order???
    grid.point_data.scalars.name = 'scalars'
    grid.dimensions = data.shape

    iso = tvtk.ImageMarchingCubes(input=grid)
    mass = tvtk.MassProperties(input=iso.output)
    return mass.volume

main()

This yields:

Ideal volume 14.1371669412
Coarse approximation 14.7969924812 Error 4.66731094565 %
VTK approximation 14.1954890878 Error 0.412544796894 %

edited Jan 18, 2012 at 14:41

answered Jan 18, 2012 at 6:42

Joe Kington

287k73 gold badges621 silver badges474 bronze badges

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

2 Comments

tylerthemiler Over a year ago

This looks awesome! I will try it out when my boss gets back (I don't even have admin privs here, so I can't install vtk :( ).

SkyLeach Over a year ago

If you get traits.trait_errors.TraitError: The 'input' trait of an ImageMarchingCubes instance is 'read only'. see my update to this answer below (code update for tvtk changes)

Charlie Martin · Accepted Answer · 2012-01-18 00:52:23Z

1

The voxels are going to be fairly simple, regular polyhedra, no? Compute the volume of each one and sum them.

answered Jan 18, 2012 at 0:52

Charlie Martin

113k27 gold badges198 silver badges268 bronze badges

1 Comment

tylerthemiler Over a year ago

The main problem for that is that the z-direction is too coarse. We do a kind of shells method to average between the slices, but this is not good enough.

SkyLeach · Accepted Answer · 2018-06-08 17:20:09Z

if you are trying to use Joe's answer above, you'll get:

`traits.trait_errors.TraitError: The 'input' trait of an ImageMarchingCubes instance is 'read only'.`

Here is the the required change and a diff showing how to fix it.

Modified Code:

import numpy as np
from tvtk.api import tvtk
from tvtk.common import configure_input


def main():
    # Generate some data with anisotropic cells...
    # x,y,and z will range from -2 to 2, but with a
    # different (20, 15, and 5 for x, y, and z) number of steps
    x, y, z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j]
    r = np.sqrt(x**2 + y**2 + z**2)

    dx, dy, dz = [np.diff(it, axis=a)[0, 0, 0] for it, a in zip(
        (x, y, z), (0, 1, 2))]

    # Your actual data is a binary (logical) array
    max_radius = 1.5
    data = (r <= max_radius).astype(np.int8)

    ideal_volume = 4.0 / 3 * max_radius**3 * np.pi
    coarse_volume = data.sum() * dx * dy * dz
    est_volume = vtk_volume(data, (dx, dy, dz), (x.min(), y.min(), z.min()))

    coarse_error = 100 * (coarse_volume - ideal_volume) / ideal_volume
    vtk_error = 100 * (est_volume - ideal_volume) / ideal_volume

    print('Ideal volume', ideal_volume)
    print('Coarse approximation', coarse_volume, 'Error', coarse_error, '%')
    print('VTK approximation', est_volume, 'Error', vtk_error, '%')


def vtk_volume(data, spacing=(1, 1, 1), origin=(0, 0, 0)):
    data[data == 0] = -1
    grid = tvtk.ImageData(spacing=spacing, origin=origin)
    grid.point_data.scalars = data.T.ravel()  # It wants fortran order???
    grid.point_data.scalars.name = 'scalars'
    grid.dimensions = data.shape

    iso = tvtk.ImageMarchingCubes()
    configure_input(iso, grid)  # <== will work
    # iso = tvtk.ImageMarchingCubes(input=grid)
    mass = tvtk.MassProperties()
    configure_input(mass, iso)
    # mass = tvtk.MassProperties(input=iso.output)
    return mass.volume


if __name__ == '__main__':
    main()

Diff of changes I made

2a3,4
> from tvtk.common import configure_input
> 
6c8
<     # x,y,and z will range from -2 to 2, but with a 
---
>     # x,y,and z will range from -2 to 2, but with a
8c10
<     x,y,z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j]
---
>     x, y, z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j]
11c13,14
<     dx, dy, dz = [np.diff(it, axis=a)[0,0,0] for it, a in zip((x,y,z),(0,1,2))]
---
>     dx, dy, dz = [np.diff(it, axis=a)[0, 0, 0] for it, a in zip(
>         (x, y, z), (0, 1, 2))]
24,26c27,30
<     print 'Ideal volume', ideal_volume
<     print 'Coarse approximation', coarse_volume, 'Error', coarse_error, '%'
<     print 'VTK approximation', est_volume, 'Error', vtk_error, '%'
---
>     print('Ideal volume', ideal_volume)
>     print('Coarse approximation', coarse_volume, 'Error', coarse_error, '%')
>     print('VTK approximation', est_volume, 'Error', vtk_error, '%')
> 
28c32
< def vtk_volume(data, spacing=(1,1,1), origin=(0,0,0)):
---
> def vtk_volume(data, spacing=(1, 1, 1), origin=(0, 0, 0)):
31c35
<     grid.point_data.scalars = data.T.ravel() # It wants fortran order???
---
>     grid.point_data.scalars = data.T.ravel()  # It wants fortran order???
35,36c39,44
<     iso = tvtk.ImageMarchingCubes(input=grid)
<     mass = tvtk.MassProperties(input=iso.output)
---
>     iso = tvtk.ImageMarchingCubes()
>     configure_input(iso, grid)  # <== will work
>     # iso = tvtk.ImageMarchingCubes(input=grid)
>     mass = tvtk.MassProperties()
>     configure_input(mass, iso)
>     # mass = tvtk.MassProperties(input=iso.output)
39c47,49
< main()
---
> 
> if __name__ == '__main__':
>     main()

Itemized list of changes:

Ran it through 2to3 so it can be run in python 3
Fixed the code for PEP8 compliance according to autopep8 (syntax, length, spacing changes)
Import configure_imput due to These changes in TVTK (Github code change)
Modify the input= kwarg in the ImageMarchingCubes constructor
Modify the input= kwargs in the MassProperties constructor
Wrap the call to main() in a direct call (to prevent execution if imported) #BestPractices

Collectives™ on Stack Overflow

Good algorithm for computing volume or surface area in python

3 Answers 3

2 Comments

1 Comment

if you are trying to use Joe's answer above, you'll get:

`traits.trait_errors.TraitError: The 'input' trait of an ImageMarchingCubes instance is 'read only'.`

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

2 Comments

1 Comment

if you are trying to use Joe's answer above, you'll get:

traits.trait_errors.TraitError: The 'input' trait of an ImageMarchingCubes instance is 'read only'.

Comments

Your Answer

Sign up or log in

Post as a guest

Related

`traits.trait_errors.TraitError: The 'input' trait of an ImageMarchingCubes instance is 'read only'.`