Ed Brignole asked about setting a volume contouring level so that the
resulting surface encloses a desired volume (cubic Angstroms). Attached
is a Python script that does this in two ways. The first way is fast
but approximate finding the level such that the number of enclosed grid
points times the voxel volume equals the requested volume. Estimating
the enclosed volume of the surface this way can give a much different
result than actually computing analytically the volume within the
surface. The second approach computes surfaces at different contour
levels and calculates analytically the volume. It uses bisection to
find the correct contour level. This approach can be slow, about 30
seconds for a 256^3 volume (sfv.mrc virus map) with ~2.5e6 surface
triangles. Most of that time (80%) is actually computing the enclosed
volume. That should be very fast but the code is checking for holes in
the surface in a rather slow way. Both methods use the step size
(subsampling) set in the volume dialog.
Tom
# -----------------------------------------------------------------------------
# Find threshold level so that a contour surface encloses a specified volume.
#
# Method: Estimates how many grid points should lie within
# surface then finds the level having that many grid points with higher values.
# This method can produce a level where the enclosed volume is very
# different from the requested volume.
#
def level_enclosing_grid_points(v, vol):
(x0,y0,z0), (x1,y1,z1) = v.xyz_bounds()
tvol = (x1-x0)*(y1-y0)*(z1-z0)
r = 1.0 - vol/tvol
s = v.matrix_value_statistics()
level = s.rank_data_value(r)
return level
# -----------------------------------------------------------------------------
# Find threshold level so that a contour surface encloses a specified volume.
#
# Method: Computes surfaces for different levels and resulting volumes
# using bisection to determine the correct level.
#
def level_by_bisection(v, vol, iterations = 10):
m = v.matrix()
l0, l1 = m.min(), m.max()
ijk_step = v.region[2]
from numpy import product
voxel_volume = product(ijk_step) * product(v.data.step)
mvol = vol / voxel_volume
return find_level(m, l0, l1, mvol, iterations)
# -----------------------------------------------------------------------------
#
def find_level(m, l0, l1, mvol, iterations):
l = 0.5*(l0+l1)
if iterations <= 1:
return l
from _contour import surface
varray, tarray = surface(m, l, cap_faces = True)
from _surface import enclosed_volume
vl, hole_count = enclosed_volume(varray, tarray)
varray = tarray = None
if vl == mvol:
return l
if vl > mvol:
l0 = l
else:
l1 = l
return find_level(m, l0, l1, mvol, iterations - 1)
# -----------------------------------------------------------------------------
#
vol = 3.227e6 # sfv.mrc, threshold level 1000.
from chimera import openModels
from VolumeViewer import volume_list
for v in volume_list():
level1 = level_enclosing_grid_points(v, vol)
level2 = level_by_bisection(v, vol)
print v.name, level1, level2
